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REMARKS 

Information Disclosure Statement 

The Examiner has requested detailed information regarding the FracCADE package 
referred to in the specification and drawings. Applicants have included herewith an excerpt from ' 
Reservoir Stimulation - Third Ed., Economides, M J. and K.G. Nolte, Wiley and Sons, 2000. 

K . " ~~ 

This excerpt describes approach used in the FracCADE package, particularly Section 6-3.2. 

The Examiner has requested information pertaining to the "well known commercially 
available simulator" described in Mr. Siebrits affidavit. Attached is a copy of SPE 80935 which 
describes the GOPHER software system. This is the system referenced in the aforementioned 
affidavit. 

The Examiner has r equested informa tion pertaining to the Lin and Keer three layer test. 
Applicants have included herewith an article from the Journal of Applied Mechanics (56 J. Appl. 
Mech.63-69) which describes this test. 



The Examiner has made several objections to the specification. Applicants have submitted 
herewith substitute pages 4, 22 and 25 to address these objections. In addition, Applicants have 
amended equation (9). 

Rejections Under 35 U.S.C. §112 

The Examiner has rejected claims 1-4, 6-10 and 12-19 under 35 U.S.C. §112, first 
paragraph as containing subject matter which was not described in the specification in such a way 
as to enable one skilled in the art to make or use the invention. 

Applicants respectfully traverse the rejection. Applicants have amended equation (9) to 

18 



Attorney Docket 56.0468 

indicate that amend the term "a/* to "a j" to indicate that it is layer dependent. 

With regard to independent claim 6 and dependent claims 4, 7-10 and 12, the Examiner 
asserts that only partially active elements are updated. Applicants note that fully active elements 
include a global mass element and therefore do not require updating. 

With regard to the rejection of claims 8 and 9, both claims have been amended to require 
that a time history of fluid volumes for pumping, fluid properties and logs are required in the first 
data set. 

For these reasons, Applicants respectfully request that the Examiner withdraw this 
rejection of claims 1-4, 6-10 and 12-19. 

Claims 3, 6-9 and 12 stand rejected under 35 U.S.C. §112, second paragraph as being 
indefinite for failing to particularly point out and distinctly claim the subject matter which 
Applicants regard as their invention. Applicants respectfully traverse the rejection. 

Claim 3 has been amended to adapt the language suggested by the Examiner in paragraph 
1 1 of the office action. 

Claims 6-9 and 12 have been amended to correct the lack of antecedent basis for the 
limitation "the recalculation of fully active elements." 

For these reasons, Applicants respectfully request that the Examiner withdraw this 
rejection of claims 3, 6-9 and 12. 

Rejections Under 35 U.S.C. §103 

Claims 1-3 and 13-19 stand rejected under 35 U.S.C. § 103(a) as being unpatentable over 
GOPHER in view of Linkov, et al. Applicants respectfully traverse the rejection. 

At the outset, Applicants note that Linkov was never coded into a simulator. Clearly, this 
is an important step in producing a system capable of modeling fractures. More importantly and 
as the Examiner may know, Linkov is not capable of modeling a formation where a hydraulic 
fracture touches or approaches an interface. This is a fundamental requirement for a hydraulic 
fracture simulator. 
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In addition and as the Examiner will note from the material supplied which details 
GOPHER, this simulator is based on equations which are incorrect for the present application. 

In summary, for reasons detailed above, it is submitted that all claims now present in the 
application are allowable. Accordingly, allowance of all claims is submitted to be in order. Such 
action is respectfully requested. 



The Commissioner is hereby authorized to charge or credit any fees to Deposit Account 
04-1579(56.0468). 



SCHLUMBERGER TECHNOLOGY CORPORATION 

1 10 Schlumberger Drive, MD-1 

Sugar Land, Texas 77478 

281.285.4524 

281.285.8569 (fax) 



Respectfully submitted, 




Stephen Schlather 
Reg. No. 45,081 



20 



Geertsma and de Klerk also extended the model 
to include fluid leakoff, following Carter's (1957) 
method. Fluid loss is incorporated by assuming that 
it has no effect on fracture shape or pressure distri- 
bution. The volume of a two-wing KGD fracture is 



(6-33) 



Performing a volume balance and solution proce- 
dure similar to that of Carter, they obtained 



L = -2£^L* , erfc(5) + 4-5-1 |, 



where 



5 = 



(6-34) 



(6-35) 



To include the effects of spurt loss S p , w w should 
be replaced by w w + (%/n)S p , which is equivalent to 
the Carter relation with w replaced by vv - 
w = kw/4. 



2S P and 



Assumptions of the PKN and KGD models 

Both the PKN and KGD models contain a number 
of assumptions that are revisited in this section. 
They assume that the fracture is planar (i.e., that it 
propagates in a particular direction, perpendicular 
to the minimum stress, as described in Chapter 3). 
They also assume that fluid flow is one-dimen- 
sional (ID) along the length of the fracture. In 
the case of the models described, they assume 
Newtonian fluids (although Perkins and Kern also 
provided solutions for power law fluids), and 
leakoff behavior is governed by a simple expres- 
sion derived from filtration theory (Eq. 6-14). The 
rock in which the fracture propagates is assumed to 
be a continuous, homogeneous, isotropic linear 
elastic solid; the fracture is considered to be of 
fixed height or completely confined in a given 
layer; and one of two assumptions is made con- 
cerning the length to height ratio of the fracture — 
i.e., height is large (KGD) or small (PKN) relative 
to length. Finally, the KGD model includes the 
assumption that tip processes dominate fracture 
propagation, whereas the PKN model neglects 
fracture mechanics altogether. 

Since these models were developed, numerous 
extensions have been made that relax these 
assumptions, the most important of which are the 
solutions for power law fluids. These two models 



are still used to design treatments and are usually 
available as options in simulators. 

Similar solutions can be derived for radial frac- 
tures (see Sidebar 6C). 

: 6C;: Radial fracture geometry models ; ; : 

: Both 1 Rerkjns and Kern (1961) and, Geertsma and de Klerk . 
; (1 969) considered radial fractures; which grow unconfihed 
I |from a point source. This model is applicable when there are ; 

ho gamers constraining height growth -or when a horizontal 
ifracture is created. ■ - r " • '"I- " '.^-M 

M "'Geertsma and de Klerk formulated the radial model using t 

the same aipuments qutlined in "Derivation of Ihe : ? F '. 
iKhristianovic 6-7). The 

Jf rapture, width \s K ./■ \f\": V, \ : ..' ■ - - • ••/ 



]/:and ihe|radial ; length Rjs ; ; : ; : ; ; 



q f (4y,^5ft) 



- VwherS||:^:-: 



e^erfc(S)V^^>:i,|,;,i 



;, : , : .-::vV ; ; 5 iJ .-15C £! V7ir,..: : 



M6C-2)!; 



(6C-3):;; 



. An explicit relationship; for, pressure can be derived by con; " 
' sjdering the solution for flow from a point; source, in which 
• pase the pressure in the fracture is a function of the expres- : 
; Si6h ln(^^, where r w is the radius of the wellboreVi • ■ Jifc /<t 
QT^ model are " 



vw: = 2:17 



f? = 0.52 



m 



:%(6jC-4):: 



it- : ^ 

W-X^- 'Si-- (6C-6) : ir 



An expression for width in the case of large' fluid loss was i 
^not provided but can be found from Eqs. 6C- 1 and 6C-6.. '■'[ 



6-3. Three-dimensional and pseudo- 
three-dimensional models 

The simple models discussed in the previous sections 
are limited because they require the engineer to spec- 
ify the fracture height or to assume that a radial frac- 
ture will develop. This is a significant limitation, 
because it is not always obvious from logs and other 
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data where or whether the fracture will be contained. 
Also, the fracture height usually varies from the well 
(where the pressure is highest) to the tip of the frac- 
ture. This limitation can be remedied by the use of 
planar 3D andf>seudo-3D(P3D) models. 

The three major types of hydraulic fracture models 
that include height growth are categorized according 
to their major assumptions. 

• General 3D models make no assumptions about the 
orientation of the fracture. Factors such as the well- 
bore orientation or perforation pattern may cause 
the fracture to initiate in a particular direction 
before turning into a final preferred orientation (per- 
pendicular to the far-field minimum in-situ stress). 
Simulators incorporating such models are computa- 
tionally intensive and generally require a specialist 
to obtain and interpret the results. They are most * 
applicable in research environments, for which they 
are used for studying details of fracture initiation 
and near-well complexities such as those discussed 
in Section 6-8, rather than overall fracture growth. 
One example of such a study was published by 
Brady et al (1993). These models are not discussed 
further in this volume. 

• ?PlahafGD v modeils are based on the assumption that 
the fracture is planar and oriented perpendicular to 
the far-field minimum in-situ stress. No attempt is 
made to account for complexities that result in devi- 
ations from this planar behavior. Simulators based 
on such models are also computationally demand- 
ing, so they are generally not used for routine 
designs. They should be used where a significant 
portion of the fracture volume is outside the zone 
where the fracture initiates or where there is more 
vertical than horizontal fluid flow. Such cases typi- 
cally arise when the stress in the layers around the 
pay zone is similar to or lower than that within the 
pay. This type of model is described in more detail 
in Section 6-3.1. 

•iPSD'mddels^attempt to capture the significant 
behavior of planar models without the computa- 
tional complexity. The two main types are referred 
to here as "lumped" and cell-based. In the lumped 
(or elliptical) models, the vertical profile of the frac- 
ture is assumed to consist of two half-ellipses joined 
at the center, as shown in Fig. 6-4. The horizontal 
length and wellbore vertical tip extensions are cal- 
culated at each time step, and the assumed shape is 
matched to these positions. These models make the 




Figure 6-4. Conceptual representation of the lumped model. 



inherent assumptions that fluid flow is along 
. streamlines from the perforations to the edge of the 
ellipse and that the streamlines have a particular 
shape, derived from simple analytical solutions. 
Cell-based models treat the fracture as a series of 
connected cells. They do not prescribe a fracture • 
shape, but generally assume plane strain (i.e., each 
cell acts independently) and do not fully couple the 
calculation of fluid flow in the vertical direction to 
the fracture geometry calculation. 

In the fixed-height models described previously, no 
consideration is given to the layers surrounding the 
fractured zone. The planar and P3D models use data 
about the properties of the surrounding zones to pre- 
dict the rate of growth into these zones. For the pur- 
pose of this chapter, planar 3D models are defined as 
those in which calculation of the full 2D fluid-flow 
field in the fracture is coupled to the 3D elastic 
response of the rock, and P3D models are defined as 
those that approximate either the coupling or the 3D 
elasticity in some manner. 

Regardless of which type of model is used to 
calculate the fracture geometry, only limited data are 
available on typical treatments to validate the model 
used. For commercial treatments, the pressure history 
during treatment is usually the only data available to 
validate the model. Even in these cases, the quality of 
the data is questionable if the bottomhole pressure 
must be inferred from the surface pressure. The bot- 
tomhole pressure is also not sufficient to uniquely 
determine the fracture geometry in the absence of 
other information, such as that derived from tiltmeters 
and microseismic data (see Sidebar 6D). If a simulator 
incorporates the correct model, it should match both 
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; 6D. Field determination of fracture geometry' 

• / Fracture geometry can be determined by using the two tech- v 
niques of rhicroseismic activity and tiltmeters; Microseisms 

: can be used to locate the fracture, thus providing estimates of 
its length 1 and height, whereas tiltmeters can provide informa- 

; . ' tion about fracture width. '] r ^'r : ^y. 1t^[A'^.. : ■. ;; ; ; 

Microseisms 

Although all modeis of hydraulic fracturing assume, that the 
rock is a continuous medium, it is well known that reservoirs i 
have natural fractures, bedding planes and other weakness .. 
features that respond as a noncontinuum. Such features have = 

,. been used to image hydraulip -f ractures using seismic tech- 7 ; 

k '; ; niques. ••• ' 1 •• " -- : ■ r ; ;j • '•: 4; ^ -• ■': ; ;i " '•' - j • • ff. 

. Hydraulic fractures induce two large changes in the reser- ; 
voir as they are created. The stress in the surrounding rocks 
: is perturbed because of fracture opening, and the pore pres- 
sure is increased as a result of leakoff. of the high-pressure. ; i 

■ t fracturing fluid. Both of these features can result in the gener- ^ 

f ation of large shear stresses on many of the weakness planes- 
near the hydraulic fracture, resulting in' small shear slippages 
called microseisms or micrqearthquakes. 

Microseisms generate seismic waves that can be detected 
by sensitive seismic receivers in nearby wells. As shown in) 
Fig. 6D-1 , both compressional waves (P-waves) and shear. 

. : i waves (S-waves) can be generated by the microseism, and ■ 
these two waves travel with different velocities; If a receiver • 
canidetect both the P- and 5-waves, the time separation can 

r be determined andlhe distance to the source inferred from 



(6D-1) 



where up and u$ are the compressional and shear velocities, 
respectively, and k and are the shear and compressional : [ 
arrival times." ( • -O- ^v:' • - : .C' - : : ; " : /*' ■ ; .'. 
v The.direction in space can be determined by using a tri- L 
axial receiver to examine the amplitude of the P-wave. The:: 
P-wave has the characteristic that its particle motion (how the ; 
rock mass vibrates) is aligned with the direction of travel of 
the wave. By obtaining the orientation of the resultant ampli- ; 
tude vector at any time, the microseism can "be traced back to 
its source: - : M \v , -/V. * v ** [ , '.. 



' With multiple seismic receivers, triangulation . techniques .. 
can be employed and greater accuracy obtained. With either : 
'approach, however, the objective is to locate the zone of i 
microseisms surrounding the hydraulic fracture and deduce 
the size and shape of the fracture from this information.- 

Downhole tiltmeters 

Width development in a hydraulic fracture results in elastic .-.,'■[ 
deformation of the formation. /This deformation can be used 
for fracture diagnostics' to provide significant information about 
fracture height and width and also about formation character- •; 

istics.. * » • • ; .' ; "• ^ / • 

As a fracture is opened, the deformation of the rock 
extends for large distances into the reservoir. Although the 
deformation is small at distances of more than a few tens of 
feet, highly sensitive tiltmeter devices can measure these 
small changes in. position. A tiltmeter does not actually mea- 
sure the ; displacement of the earth, but rather the curvature of 
the displacement, and itis capable of measuring up to nano- 
radian resolution (a nanoradian is the angle induced by 
stretching a line from New York to Los Angeles and raising ..• 
the New.York side by the diameter of; a pencil). Tiltmeters • 
have long been used for surface diagnostics of earth move- 

. ment, but the application of a string of downhole tiltmeters 
provides highly sensitive fracture data. . • • 

■ Figure 6D-2 shows a schematic of the tilt response of the \ 
formation measured in a well offset to the fracture treatment. ' 
The characteristic S-shaped curve is typical of tilt, as opposed ;. 

. to strain; and can be simply explained. Straight across from ; 
the fracture, the rock is -pushed away, but is not tilted on the 
geometric axis of the fracture, and there is zero tilt. Above the 
fracture, the earth experiences curvature that is defined as o.v •. 

, negative for this example. The curvature i reaches a maximum 
at a well-defined point and then decreases to zero as the dis- i 
tance.f rom the fracture increases: The bottom is identical to 
the top/except that the curvature has the opposite direction : 
and opposite sign. /■;"■;:.!■ -^^-y 

Two aspects of this distribution are important for diagnos- 
tics.: First, the locations of the maximum tilt values are a Junc- 
tion of the height h { of the fracture relative tq, the distance : of . 
away. Thus, fracture height can be- quickly estimated. Second, ; 
the amplitude of the tilt is a function of the width of the frac : 
ture, so; the width during fracturing, and possibly the final : 
propped width, can be estimated as well. 
; : Branagah etal: (1996) provided an example of the applica- 

[ tion of tiltmeters to the;. calculation of hydraulic fracture geometry. 
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Figure 6D-1. Microseismic traces at the receiver resulting 
from shear slippage. y ^ : y ; v ' v ' i ' p' r ' -t y'.- 
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Figure 6D-2. Tiltmeter response to hydraulic fracture 
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treating pressure and fracture geometry. These issues 
are addressed in Section 6-12 and Chapter 9. 



6-3.1. Planar three-dimensional models 

A planar fracture is a narrow channel of variable 
width through which fluid flows. The fracture geome- 
try is defined by its width and the shape of its periph- 
ery (i.e., height at any distance from the well and 
length). Both the width at any point and the overall 
shape vary with time. They depend on the pressure 
distribution, which itself is determined by the pressure 
gradients caused by the fluid flow within the fracture. 
Because the relation between pressure gradient and 
flow rate is highly sensitive to the fracture width 
(Eq. 6-9), the geometry and fluid flow are tightly cou- 
pled.. Although the mechanics of these processes are 
described separately in this section, the complexity of 
solving any fracture model lies in the close coupling 
between the different processes. Three separate prob- 
lems are considered: 

• width profile in a fracture of known shape and pres- 
sure distribution 

• shape of the fracture 

• flow of fluid in a fracture of known shape and width 
(i.e., known geometry). 

Hirth and Lothe (1968) and Bui (1977) showed 
how the pressure and width in a fracture may be 
related. Basically, the width at any point (x,y) is 
determined by an integral of the net pressure over 
the entire fracture, expressed as 

= JJ /(* -*',y- y')(p(*', y) - a(x'. y'))dx'dy\ 

s 

(6-36) 

where a is the stress. 

The details of the elastic influence function / in 
Eq. 6-36 are beyond the scope of this volume. Use- 
able forms of Eq. 6-36 can be derived generally only 
for homogeneous linear elastic materials (see Side- 
bar 6E). In fracturing applications, the rock is usually 
also assumed to be isotropic. 

The shape of the fracture evolves with time. In 
essence, the boundary (i.e., the vertical and horizontal 
tips) moves outward as the fluid provides sufficient 
energy to fracture the rock at the boundary. More com- 
plex tip behavior is discussed subsequently, but in this 



section it is assumed that this process is described by 
linear elastic fracture mechanics (LEFM). If the LEFM 
failure criterion is exceeded at any point on the fracture 
periphery, the fracture will extend until the criterion 
is again met. For simple shapes and pressure distribu- 
tions, such as ellipses under constant pressure, the cri- 
terion can be specified analytically, similar to Eq. 6-3. 
For more complex shapes and pressure distributions, 
analytical solutions are not available. In these cases, 
it can be shown that a relatively simple criterion can be 
written in terms of the width near the tip and the criti- 
cal stress intensity factor or fracture toughness K fc , 
which is introduced in Chapter 3: 



w(x) = - 



(6-37) 



where x is the distance measured from the tip. Rela- 
tions between fracture mechanics parameters such as 
the specific surface energy (used in Eq. 6-3) and the 
fracture toughness are provided in Chapter 3. 

The fluid flow is described by equations for conser- 
vation of mass (a general form of Eq. 6-21, including 
the density p and expressed in terms of velocity u): 



dx 



dy 



+— (pvv) + 2p W ,=0, (6-38) 



which can be written as a vector equation; 

V(pvvw) + ^-(pvv) + 2pw^ =0, 
at 



(6-39) 



and the conservation of momentum (a general form of 
Eq. 6-9) is 



P^ = -V/>-[V-T] + pg, 



(6-40) 



where x is the shear stress and g is the acceleration of 
gravity. 

The first two terms in Eq. 6-38 relate to the spatial 
change of the mass-flow vector, and the second two 
terms represent the storage resulting from width 
increases and leakoff, respectively. Equation 6-40 is 
a vector equation. The term on the left-hand side is the 
rate of change of momentum, and the terms on the 
right-hand side are the pressure, viscous and gravita- 
tional forces, respectively. It simply states that a small 
element of fluid accelerates because of the forces act- 
ing on it. This equation can be expanded and then sim- 
plified for the geometries of interest in hydraulic frac- 
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6E. Lateral coupling in pseudo-three- : - 
:-;vV ; " dimensional models^,: • i" : -!"-:.: 5 ^' f^/- ' 

Assume that a. fracture has a fixed height and' that it consists : f . :i 
; of a number of elements each of constant width over the height 

(i.e., a KG D fracture). Let the grid points be represented by :: ;i 
. points X/ in the center of the elements with comers- (x^y^tj, [ . . 

,' Start ield ( 1983) developed a boundary element solution tech : . ;: ' 
nique called the displacement discontinuity method. They -.'/J 
showed that the pressure at any point is given by. ' ; ;. : ■'' : 



Yt.i = Yt,k 











— wv 


(x Cl i.yc,i) _ 








— ^(x Cl k,y C( k) 










— w 



I 

I 

I 

x l.i 



yb,i = yb.k 



i 



x J,k *r,k 



ice function of 



where ^jc: is an influence function of the form : 
where the influence 'function / is defined as 



:(6E-1) ; : : 



(6E-2)i ;. 











tiy.cj-y,,) 2 




[(xj+x-tfi 




jV2 ; :^ 


Wcy-Xrj,) 2 :* 


)(y e ,+y»J 

vfycy+'y^) 2 : 





•(^/;-^,v)(yci+yM.) 



;i6E-3): 



i; lb accurately solve Eq: 6E-1 requires a large number of : - r 
elements. Also, it is-difficulttd extend directly to other shajpes • 

• such -as ellipses ;br for nonconste^ 
these problems; the equation is modi^ as follows. The ; ; : 
wiolh at any point can be written as V : : - :; ;\ ; ; | : 

• ; where Aivw.is defined £s }fi^ : s^;,j;;i s-M^-: '■■ V:vi-- 



-Figure 6E-VG^ 



! i J ^Equation %E-1 can then be written .asyi.ii- - ; ': i -f?:. : r 
^ : ' ! ^whe're|.::: : :; ; •; ; ;'^Ji r t:: -^l^'i ^5:CV. -^fr; r ; j: : V' ; :'v ',. . 

, The term ^3c)X^* thus represents; the - pressure induced 1 i; 

:by a fracture of constant width^ 
: : length, this pressure would be exact if calculated using the f :^ :} 

plane strain solution - TTjeiitemij^ w) can therefore be \ r -W\ 

obtained as the sum of the plane strain solution and the effect 

of two semi-infinite fractures of vyii- ^ attached at the tip of i 
:. : : k each fracture wing. '{ MM :[ 2 

. From: Eq • 6E-2, the influence fuhctionadecrease: with ^ % ; 

distance from an element The advantages of the form of ' ; t'pi 
: ; : "Eq.; 6E^8 are that the correct ■ 
$ rWeiht where the widths are almost the same and that the :| ' ; I : if J 
; ; self-correction is exactly zero by definition. The number of • 
; • ; elements riequired ; to obtain an accurate solution is significantly ; 
; reduced,^ and variable heights ;ahd; other shapes are easily - V • 

introduced. Lateral coupling is relatively easy to introduce to^:; 1 S 
1 : the explicit solution method because the pressure .correction ' ; 
v is simp|y. added before the:^ are calculated : • '/.:f : 



turing (see Sidebar 6F). For a particular component, 
such as the x component, Eq. 6-40 can be written as 

dx 



Du,. 



Dt 



dp 
dx 



. dr., 



dy dz 



+ P&. (6-41) 



A constitutive law relating the stresses x to the flow 
rate is required to complete the description of fluid 
flow. In the case of steady flow in a narrow channel 
such as a fracture, the full details of the constitutive 
law are not required, because the narrow fracture 
width results in the complete dominance of some 
stress terms. The only terms of interest are the shear 
stresses induced by velocity gradients across the frac- 
ture. In addition, use is made of the lubrication 
approximation, so flow perpendicular to the fracture 
wall (the z direction) is neglected. With these assump- 



tions, the equations for the stress in a Newtonian fluid 
reduce to 



t - = = Ait 

x » = x - = -is 



and Eq. 6-41 can be written as 



Du. 



dp 



2„ \ 



d 2 u 



(6-42) 



(6-43) 



For the special case of a narrow channel (Poiseuille 
flow), where velocity gradients parallel to the flow are 
small and there is no flow perpendicular to the chan- 
nel, the time-dependent term simplifies to a partial 
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6F. Momentum conservation equation for : ; 
hydraulic fracturing 

Equation 6-40 is a vector equation, for which a component 
: : can be written' 'as ' : : "Cv ■ ' " : . •.' ; . , . -m,- 1 



Dt 



(6F-1) 



where t/ is the velocity, cms the gravitational acceleration, and 
Vis x, y or z. The term on the left side of Eq. 6F-1 is termed 
the substantial derivative, which is the rate of change seen by 
an observer moving with the^r can be related to 

the usual partial derivative (i.e., the rate of change seen by a 
^stationary observer) as^; ! ; - "-oV;^^ 'f If^i ;.';'" ' : 



' ■ d a a a ~ i .. a .; 

:■ v.. - . . .; =; — + +u — : - : 

. ; ; - , KDt . . ;dt ; ' ax : ? ; dy -J* :'-dz.,-- ■ 
Thus, Eq. 6F-2 can be expanded to . .: r ; - 1; ; 



(6F-2) 



dx, 



ax ay . dz 



(6F-3) 



; : ;V. This completely general equation can be simplified for a. 
; narrow channel in an impermeable medium. Leakoff does not 

occur in this case, so components in the z direction can be 
; jneglected. In addition/ the flow is assumed to be steady state, 

so time derivatives can be ign^ 
• plifies to : i .'. .. " ]\'- ' ; •/ ; ; ^ V''=C \ :: •' '• '.' . ; " ■ 

; ^i^"f ; ; : • au/ - aii,.0 I Milf ■ \-r ■ " v - ;; - : v ' 



ax,. : 



^ ax; ; ' c)y : ^Z; 



(6F-4) ; 



for 7.= 1 or 2. Even for a permeable medium, Eq. 6F-4 is 
used. In this case, leakoff is treated as a sink term and 

, . included in the mass balance, but it is assumed not to affect ; 

; the equations relating pressure, stress and fluid velocity; 

: Newtonian fluids : ; / : - . y v '/' V' 

!|' To make Eq. 6F-4 useful, the stress components must be 
determined, which is done by assuming a model of fluid , 
; behavior: For example, a Newtonian fluid is a'model with one 
V parameter, the viscosity ji. The stress components^ ^ ; . 

'•• - ; k ■ : , ■ ■ .■: : -:du-- V z 2 ■■ yly v ' \ •'•',/" ' '• \ ' . 



(6F-5); 



The first: three components of Eq. 6F-5 are the normal 
stresses, and the last three are the shear stresses. The last-, 
term of the normal components is zero for incompressible 
fluids. In the case of 1 D flow between parallel plates, without; 
leakoff, two of the velocity components are identically zero. 
In addition, conservation bf mass implies that the third com- • 
ponent cannot vary with position. Hence, all three normal . ' 
components are identically zero. The equations thus reduce : ' . 
to those' for shear flow. Although these assumptions are not ; • 
strictly true in general, they are used for the flow calculations 
in hydraulic fracture modeling. It can also be shown that for • 
a narrow channel, the velocity- gradients perpendicular to the ■ 
walls (the z direction) are much greater than those in the par- 
allel directions. Finally, therefore, the stress components for J 
a Newtonian fluid jn a hydraulic fracture can be written as ; :• 



Substituting Eq. 6F-6 into Eq. 6F-4 obtains : : 

i . ' ( aJ, i : ': . au, ¥^ dp •• Wu J'-fl 



P \M, 



dp :d°u r 
ay J By - H az 2 



■For 1 D flow along the fracture length, as typically assumed in 
P3D models, Eq. 6F-7 can be simplified to .4 



•az 2 / 



:Yap : 
p. ax " 



(6F-8) 



^Assuming zero slip (i.e., zero velocity at the fracture wall), the 
solution to Eq. 6F-8 is : i : } •. . ■ ^}[' ' '^- : y :K 



: (6F9) 



^Integrating to obtain the average velpcity across the channel, 

: .'' : ' :V ;: L . (6F-10); 



- _ -w r dp 

~\ ,l2ji-. dx' 



The flow rate; per-unit height is obtained by multiplying the ' 
; average velocity by the width w. , ' : , : :> . ; ; ;^ ; - •* - • ; , : : \ : :J 
In the case. of 2D flow; the left-hand sides of Eq. 6F-7 are; ^ 
zero if inertia may be neglected. In this case for the y direct 
tion, an equation can belformed similar to; Eq. 6F10, except ; 
that it includes a gravitational term. ;; ;;;- ■ • .'; ';'4;_ 
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derivative of velocity with respect to time. It is usually 
assumed that the flow is steady state, which finally 
obtains 

'3V 



dz 2 



(6-44) 



and a similar equation for the y component. 

Equations 6-36 through 6-44 summarize the planar 
3D model for Newtonian fluids. Similar results can be 
obtained for non-Newtonian fluids (see Sidebar 6G). 
These equations are generally not amenable to ana- 
lytic solutions, but require a numerical simulation. 
In addition, although it is relatively straightforward 
to write the conceptual equations, efficient and robust 
numerical solutions are difficult to obtain. The pri- 
mary reasons for this difficulty are the extremely 
close coupling of the different parts of the solution 
(e.g., fluid flow and solid deformation), the nonlinear 
relation between width and pressure, and the com- 
plexity of a moving-boundary problem. 

The first numerical implementation of a planar 
model was reported by Clifton and Abou-Sayed 
(1979). In essence, their approach was to define 

6G. Momentum balance and constitutive equation : = 
for non-Newtonian fluids ' ' : --.A-"f 

, : ; The definition of a Newtonian fluid is the one-parameter rela- \ 
tion between stress and velocity {Eq. 6G-5). In tensor notation, -A, > 
this.can be written as ,' ^ . 

- : f : • =-nA^: Y - : • (6G-1); 

• where A is the rate of deformation tensor, with components I f 



_ dUi ^ 



du 



(6G-2): 



The viscosity may be a function of pressure and tern ] 
ture or other variables, including the history of the fluid, but not i 
of A. For non-Newtonian fluids, ah equation similar to Eq. 6G-1 
may be written: a, "A : " v = r'AA'-;K- ■■ 



■(6G-3);: 



where |w is a function of a; For flows of the type of interest 
in fracturing, it can be shown that u.« may depend only on A 
through a relation of the form <. -Ay. -a '. j' ; 

'.'■,./ • ; [A.AA ' . (6G-4) : 

where, h is the second tensor invariant: , 

':^:£-:\a V = XE a * a V- .■ a:- =' " 

For example, for a power law fluid, the function ji B is 



a small fracture, initiated at the perforations, divide it 
into a number of equal elements (typically 1 6 squares) 
and then begin solution of the equations. As the 
boundary extends, the elements distort to fit the new 
shape. One difficulty with such a solution is that the 
elements can develop large aspect ratios and very 
small angles, as shown in Fig. 6-5. The numerical 
schemes typically used to solve the equations do not 
usually perform well with such shapes. 

A different formulation was described by Barree 
(1983), and numerous field applications have been 
reported (e.g., Barree, 1991). It neatly avoids the prob- 
lem of grid distortion by dividing the layered reservoir 
into a grid of equal-size rectangular elements, which 
are defined over the entire region that the fracture may 
cover. In this case, the grid does not move. Instead, as 
the failure criterion is exceeded, the elements ahead of 
the failed tip are opened to flow. and become part of 
the fracture, as shown in Fig. 6-6. Two limitations of 
this approach are that 

• the number of elements in .the simulation increases 
as the simulation proceeds, so that the initial num- 
ber may be small, resulting in inaccuracy 



and for a Bingham plastic 



(6G-7); 



The commonly used consistency index K is dependent on ,, : 
the flow geometry and is related to a basic fluid property, the 
generalized consistency index K (Eq. 6G-6)- For parallel !• 
;plates;(i.e:, in a slot), 1 which can represent a fracture, the relav 
■ - tionship is A; ) ■ "A) _ jy '"--Vy ' -'Y ■' ] { Y' '.V; ■;: .■ . /' 



For a pipe it is 



(6G-8) 



(6G-9) 



The maximum difference between the two expressions is less 
than 4% for all values of n. For ID flow of a power law fluid ' 
between parallel ^plates, the average fluid velocity is given by 



n f w 
1 + 2nt,2 



(6G-10) 



For the special case of the power law exponent n = 1 , this 
reverts to the equation for a Newtonian fluid, with K replaced 
by the viscosity. Table 6G-1 summarizes useful information for 
the laminar flow of both Newtonian and power law.fluids : 
under different geometries. However, the expressions for 
pressure drop are not generally applicable for drag-reducing y 
fluids such as those used in hydraulic fracturing. : 
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6G. Momentum balance and constitutive equation for non-Newtonian fluids (continued) 



Table 6G-1. Summarized expressions for laminar flow of Newtonian and power law fluids. 



Fluid Type ; Pipe, 



Parallel Plates 



Ellipse (zero 
eccentricity) 



■Reynold's number (Afe.) 



Newtonian 



Power law 



2puw 



7tpi7lV ; 

: : :-:"2^" 1 - ; ;^" : 

• 2"/C 



; Hydrau|ic diameter (D H ) 
Friction factor .; 



16/A/ fle • : 



2w "> r , 

24/AW< 



; Velocity distribution 



Newtonian 



Power ; law 



?SSf 



^Pressure drop {&pllox:dpldx) Newtonian ;: 



Power law 



128u,<? : 

itD* r,v 



;;:/ 4n + 2 yj : 2g"/C.: 



; 64fiQ : : *0 

See Eq. 6-57 ' 



y 




























X 



Figure 6-5. Planar 3D fracture divided into elements that 
were initially square, 

• the general size of the fracture must be estimated in 
advance of the simulation to ensure that a "reason- 
able" number of elements is used. 

In addition, this particular implementation has two 
simplifying assumptions, that a simplified method is 
used for representing modulus contrasts and a tensile 
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Figure 6-6. Fixed-grid solution showing elements open to 
advance the fracture. 

strength criterion is used for fracture extension, rather 
than a fracture mechanics effect. The failure criterion 
is used to compare the stress at the center of all 
boundary elements with the material tensile strength. 
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If the strength is exceeded, then the element is 
assumed to open. However, the fracture-induced stress 
in the material near the tip of a fracture varies with the 
square root of the distance from the tip. Hence, the 
failure criterion is grid-resolution dependent. 



6-3.2. Cell-based pseudo-three- 
dimensional models 



In cell-based models, the fracture length is divided 
into a number of discrete cells. This is directly analo- 
a " (L-gous to the planar models, except that only one direc- 
\$fi a(\(^ ti° n * s discrete instead of two. Fluid flow is assumed 
0-y\ to be essentially horizontal along the length of the 
fracture, and the solid mechanics is typically simpli- 
fied by assuming plane strain at any cross section. As 
in the PKN model, these assumptions make these 
models suitable primarily for reasonably contained 
fractures, which are long relative to their height. 

These two assumptions allow separating the solid 
and fracture mechanics solution from the fluid flow 
as follows. Plane strain implies that each cross section 
acts independently of any other. In addition, the 
assumption of ID fluid flow implies that the pressure 
in the cross section is always 



P = P cP + P8y> 



(6-45) 



where p cp is the pressure along a horizontal line 
through the center of the perforations and y is the ver- 
tical distance from the center of the perforations. 
Equation 6-45 is valid only if vertical fracture exten- 
sion is sufficiently slow that the pressure gradient 
resulting from vertical flow can be neglected. This 
assumption that the vertical tips of the fracture are 
approximately stationary at all times is called the 
equilibrium-height assumption. 

• Solid mechanics solution 

With the equilibrium-height assumption, the solid 
mechanics solution simplifies to the determination 
of the fracture cross-sectional shape as a function 
of the net pressure, or p cp . Simonson et al. (1978) 
derived this solution for a symmetric three-layer 
case. Fung et al (1987) derived a more general 
solution for nonsymmetric multilayer cases. 
Following Fung et al the stress intensity factors 
at the top and bottom tips K hi and K\u respectively, 
can be written in terms of the pressure at the center 
of the perforations p cp and the closure stresses in 
the layers o, as 



fatf M '{2 { k, J 



(6-46) 



K u = 



vT 



2 



h f -2h { 



(6-47) 

where p r is the fluid density, h cp is the height at the 
center of the perforations, and hi is the height from 
the bottom tip of the fracture to the top of the /th 
layer, as shown in Fig. 6-7. 

This set of nonlinear equations can be solved by 
iteration. Assuming that the solution (two vertical 
tip positions plus the pressure) at one value of p cp is 
known, a height increment is assumed. The incre- 
mental height growth in the two vertical directions 
is then calculated such that Eqs. 6-46 and 6-47 are 
both satisfied, and p cp to obtain these positions is 
calculated. Finally, the width profile associated 
with this solution can be obtained as 




Pc P + ?fg{K P - y) " <y, )M h f~ y ) 



n-\ 

I (<*,♦, 
7 = 1 



(h, - )»)cosh" 



I y - h t I h { \y-h. 



■J y(h,-y) cos' 



hj-lh, 



(6-48) 



where y is the elevation measured from the bot- 
tom tip of the fracture. 
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Figure 6-7. Definition of variables for the fracture contain- 
ment problem. 




Consider, for example, the symmetric three-layer 
case shown in Fig. 6-8. If the gravitational compo- 
nent is neglected, so that the problem is symmetric, 
then the penetrations into the two barriers are 
equal. In this case, Eq. 6-46 can be simplified sig- 
nificantly and written as (Simonson et aL, 1978) 



K lu - K Il - 



~1 



p -a 

r C p ^ pay 



2Aa 



-cos 



71 























J 





(6-49) 

where Aa is the difference in stress between the 
central layer (pay zone) and the surrounding layers, 
and hpay and a w . are the thickness and stress of the 
pay zone, respectively. Figure 6-9 shows fracture 
height as a function of net pressure, as calculated 
by Eq. 6-49. 

Although Eq. 6-49 is for a special case, it shows 
two interesting practical results. First, penetration 
into the barrier layers occurs at a critical net pressure: 



2Kl 



(6-50) 



For example, if K lc is 2000 psi/in. 1/2 and h f \s 20 ft 
[240 in.], the critical net pressure for breakthrough 



1.5 



1.0 



0.5 



K*c= 1000psi/in. 1/2 / 
02-0! = 1000psi ^ 
_ o 2 - Oi = 500 psi / . 
c 2 -Oi = 250 psi/~" 
I I 1 




Figure 6-8. Simple three-layer height growth problem. 



Figure 6-9. Fracture height versus net pressure for sym- 
metric barriers. h 5 = penetration into the boundary layer. 
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is only about 100 psi. Second, the net pressure can- 
not reach the stress contrast because this would 
result in infinite fracture height. 

In a typical cell-based simulator, a table of these 
solutions is calculated prior to simulating the frac- 
ture evolution rather than at each time step of the 
calculation, and the relations among width, pres- 
sure and height are used to greatly speed up the 
solution of the fluid flow equations (conservation * 
of mass and momentum). 

Fluid mechanics solution 

One of the major differences between planar 3D 
and P3D models is the fluid flow calculation. The 
fluid flow model in most P3D models is the same 
as that introduced by Nordgren (1972) (i.e., a ID 
version of the model described for the planar 3D 
model). In this model, both vertical flow and the 
variation of horizontal velocity as a function of 
vertical position are neglected. This results in the 
inability of typical P3D models to represent several 
aspects of behavior, namely (Smith and Klein, 1995) 

- effect of variations in width in the vertical direc- 
tion on fluid velocity 

- local dehydration, which is approximated as 
simultaneous dehydration over the entire height 
of the fracture 

- fluid loss after tip screenouts (TSOs), when fluid 
flow through the proppant pack is ignored 

- proppant settling resulting from convection or 
gravity currents. 

The average velocity and width are used (width 
is replaced by cross-sectional area divided by 
height) to simplify the conservation of mass 
(Eq. 6-38 for an incompressible fluid) to 



dAu dA ^ / , \ 
— + — = -2V «i, , 



(6-51) 



where u is the average cross-sectional velocity and 
u L and h L are the leakoff rate (Eq. 6-14) and height 
in each layer. Similarly, the conservation of 
momentum simplifies to 

dx dz 



(6-52) 



For a power law fluid with properties n and K, 



[dz 



(6-53) 



Solving Eq. 6-52 with Eq. 6-53 with the no-slip 
boundary condition (i.e., zero velocity at the frac- 
ture wall), the average velocity across the channel is 



u x = -sgn 



fdp} 

{dx, 



dp 



dx 



K 



l+2n 



(6-54) 



where sgn represents the sign of the quantity. 

In the special case of a Newtonian fluid, n - 1 
and |i = K, and Eq. 6-54 becomes 

_ w 2 dp 

To obtain the total flow rate across the height of 
the cross section, and hence an average velocity for 
substitution in Eq. 6-51, Eq. 6-54 is integrated from 
the bottom to the top tip of the cross section: 



q = jw(y)u x (y)dy. 



(6-56) 



The average velocity is thus determined as 





( 


dp 


H 















u = — - -sgn 
A 6 



where the channel function O is 



(6-57) 



(h \ 



1 f ]+2n 

-Jw(y) « dy. (6-58) 

7 /., 



2 + 4n h 



Relations for the PKN model with power law 
fluids can be derived following this approach (see 
Nolte, 1979, 1991). 

- Laminar and turbulent flow 

When fluid flows between parallel plates at a low 
rate without leakoff, any fluid element remains 
a fixed distance from the wall of the channel, 
except in a small entrance region. This is known 
as laminar flow. By contrast, in turbulent flow, 
eddies occur, and fluid is continually mixed. This 
mixing results in added friction and different flow 
behavior. The Reynold's number N Re (defined in 
Table 6G-1) enables determining whether laminar 
or turbulent flow will occur. If N Re exceeds 2100, 
flow will be turbulent. Inside the fracture, N Re is 
typically well below this value, except for partic- 
ularly thin fluids, such as acid. 
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- Rheology of fracturing fluids 

Fracturing fluids are generally treated as power 
law fluids, and because they are shear thinning 
(i.e., viscosity decreases with increasing shear 
rate), n is usually less than 1 . The effective para- 
meters of the power law model K' and n are 
typically derived from laboratory measurements 
(see Chapter 8) over a range of shear rates. For 
shear-thinning fluids, the apparent viscosity 
(derived from K' and n) decreases as shear rate 
increases, and the viscosity would be infinite at 
zero shear rate. In reality, limiting low- and high- 
shear viscosities occur and must be considered. 

Fracturing fluid properties change with time 
and temperature. Typically, exposure to high 
temperatures reduces fluid viscosity. However, 
crosslinkers may cause initial viscosity increases 
prior to the degradation. The effects of tempera- 
ture and time are included in numerical hydraulic 
fracture simulators, typically by means of tables 
of K' and n versus time at a series of tempera- 
tures, which are similar to those in service com- 
pany handbooks. 

• Numerical solution of the model 

The three basic solutions described for height- 
growth mechanics (pressure-width-height relation), 
conservation of mass and conservation of momen- 
tum (velocity-pressure relation) are coupled and 
solved simultaneously. There are several methods 
by which the coupled equations may be solved, two 
of which are introduced here. Either a fixed or 
moving mesh may be used for the two methods, 
as described previously for planar 3D models. In 
this section, the explicit finite-difference method 
is introduced with a grid that moves with the fluid 
and an implicit method is described. In each case, 
prior to starting the simulation of the fracture evo- 
lution, a table of the pressure-height- width relation 
(from the equilibrium-height solution) is calculated 
as described for "Solid mechanics solution" in 
Section 6.3-2. 

For the explicit finite-difference method, the fluid 
in the fracture at any time is divided into n ele- 
ments, each with a cross-sectional area A t and 
bounded by two vertical surfaces at x, and x i+] , 
moving at velocities u t and u i+] , respectively, as 
shown in Fig. 6-10. (The grid is numbered such 
that i = 1 represents the tip to facilitate the addition 
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Figure 6-10. Fracture divided into elements with positions 
and velocities defined at grid points. 



of new elements at the well, as necessary.) Mass- 
conservation Eq. 6-5 1 can be rewritten as 

with the derivatives replaced by central finite- 
difference approximations to obtain 

A ' r =A; "t VL+ a; ^ - ) £ • (6 - 60) 

where V L represents the volume leaked off over the 
element of length Ax in time step Ar. The velocities 
are calculated at the grid points, and the area is 
assumed constant in each element. The cross-sec- 
tional area can thus be updated from the values of 
the velocities and areas at the previous time step. 
Once this has been done, the pressure at each cross 
section can be obtained from the solid mechanics 
solution by looking up the pressure in the precalcu- 
lated pressure-height-width relation table for the 
corresponding area A. Pressure gradients can then 
be calculated using the approximation 

*P - P»-Pi (6-61) 
d*i (*/-i-*w)/2 

and new velocities obtained using Eq. 6-57. Once 
all the velocities are known at a given time, the 
positions of the grid points are updated using 

x,.(r + Ar) = x.(t) + u.(t + A/)Af . (6-62) 

This method is known as a Lagrangian method 
because the grid coordinates move with the fluid. 
Leakoff causes each element to shrink and possibly 
even disappear as it penetrates farther into the frac- 
ture, limiting the usefulness of this method for 
modeling hydraulic fracturing treatments. In addi- 
tion, new elements must continually be added at 
the wellbore. This makes it difficult to control how 
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many elements are present at any time or the sizes 
of the elements. Another approach is to introduce a 
fixed grid, as discussed for planar 3D models. This 
has the advantage that the number of elements in 
the simulation is relatively small near the beginning 
of the simulation when less accuracy is required 
and increases as the simulation progresses. Yet 
another approach is to introduce a moving mesh 
in which the grid points move at some reasonable 
velocity, for example, such that the fracture is 
always divided into a fixed number of equal-size 
elements (i.e., using a stretching coordinate system; 
see Sidebar 6H). 

One of the primary limitations of explicit finite- 
difference methods, such as those introduced in the 



preceding text, is that the time step used in the cal- 
culation may not exceed a critical value to ensure 
stability. Because only quantities from the previous 
step are used in moving forward, numerical errors 
can grow larger from step to step if the time step is 
too large. In the development of a general hydraulic 
fracturing simulator using such differencing 
schemes, the time step must be chosen carefully to 
avoid stability problems and yet minimize the com- 
putation time. A simple stability analysis is in 
Sidebar 6H. 

It has been found that in cases of high leakoff 
or large widths (such as TSO designs), the critical 
time step for stability may be too small for efficient 
solution of the system, limiting the utility of the 



6H. Stretching coordinate system and 
stability analysis 

Stretching coordinate system 

One way to simplify grid point bookkeeping is to use a 
stretching coordinate system. If : v ; ' ■ : , f : 



X = 



(6H.1); 



then X will always remain bounded between 0 and 1 while x ; : 
varies between 0 and /.(/). Placing a grid on Xwill fully coyer \\ 
the fracture regardless of the growth: characteristics. However, 
although ^he gridding is; simplified, the complexity of the differ-: 
ential equation is increased. The derivatives are found as ! :| 



df : 



= i dt L di dX 



Equation 6-59 becomes, ';^My, 
; f . : dt L dt dX LdX • 



(6H-2); 



(6H:3) 



(6H-4) 



and the other equations, of the system are similarly trans^ . 
formed^ : : ■ ■ .. : ■■■}/;' ■ ■■ ':;.': ?":-:; •• . ; . - ; >:V-v\ •. : .-• 

Stability analysis ' 

A full stability analysis for a nonlinear system is difficult, but 
an approximate time-step limitation can be found as follows.' 
Assume that the pressure gradient can be written as 



(6H-5) :. 



In the case of the PKN r model, where the fracture height h, js 
fixed, C p = f>ht, where p is defined by 



(6H-6) 



Substituting Eq. 6H-5 into Eq. 6-59 and applying the chain, rule, " 

(6H-7): 



dt : . n 



•dx 



where absolute values must be assumed for all quantities, ; 
' : because an error.analysis is.being performed, and 0 is 
defined as- v ; ' ; v- •-• • . - '.' 0 " • '.. ' : - v ; . 



; 0 = 



the highest order term in Eq. 6H-7 is 



(6H-8); 



(6H-9) 



•If the derivative is expanded using a central difference 
approximation, the term in A becomes - 



:p-n(Ax) dx 



(6H-10) 



l To investigate the effect of an error introduced into A A is • 
replaced by A(1 .+ e), which can be approximated (for small 



[ -2D dA 
n(Ax) 2 . dx - 



• (6H-11):; 



: If a time step is taken (discretizing Eq. 6H-7 similar to Eqi 6-59), 
: then the error e grows to. ■ \ '~ ... 



-2D dA 



: n(Ax) 2 dx • 

m< 

M 



^^^(^i 3 ^ ^^) - : (6H-12) i 



For this error, to reduce in magnitude, it must be smaller than 
Ae,, which can occur only if : • . 



Af < - 



2AC P C V ■, .,:. ■ : .: 
where the viscosity leakoff control coefficient C v is 

:.:>. --\ ,; / ■ . dx \. . . • 



(6H-13) 



(6H-14) 
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explicit scheme. An implicit finite-difference 
scheme, with no time-step limitation, may elimi- 
nate this limitation. In essence, implicit and explicit 
methods can be distinguished by the fact that 
explicit methods solve for quantities at one time 
step on the basis of only values at the previous time 
steps, whereas implicit methods include the use of 
quantities at the current time. This implies that a set 
of equations is set up and solved, because the quan- 
tities at the current time step must all be found 
simultaneously. For linear problems, a set of linear 
equations results, and these are easily solved by 
standard methods such as gaussian elimination. For 
the ID flow problem, the implicit finite-difference 
formulation yields a tridiagonal system of equa- 
tions (i.e., a sparse matrix with only three diagonals 
filled with nonzeros). Highly efficient solution 
techniques are available to solve such systems . 
(e.g., Carnahan etal, 1969). For nonlinear prob- 
lems, however, such methods can be complex and 
are not always much more efficient than explicit 
methods. Iteration is frequently required, because 
a nonlinear system is linearized. If the linearization 
approximation is inaccurate, it must be corrected 
and resolved. 

Another method without the time-step limitation, 
and which avoids forming a system of equations, 
is a method using integrated or analytical elements. 
A similar method to that described in the following 
was the basis of the commercial time-sharing 
method made available by Amoco between 1981 
and 1983 (Nolte, 1982, 1988a). Consider once 
again the basic equations of the PKN model with 
x = <|> at the tip: 



w = - 



2p,...h 



net / 



E' 



dp _ 64#|i. 



nh f w* 



(6-63) 



(6-64) 



Substituting Eq. 6-64 for p into Eq. 6-63 obtains 
E' dw 64q\i 



2h f dx 



Kh f w 3 ' 



(6-65) 



Detailed numerical simulations have shown that 
the velocity varies much more slowly than the flow 
rate q because the reduction in width toward the tip 
partially compensates for fluid leakoff and storage 
in the fracture. Instead of the Perkins and Kern 



(1961) assumption that q is constant (Eq. 6-10), 
replacing q by nuhfw/4 allows writing Eq. 6-65 as 



dw _ 32(1/1, 



dx 



E'w 2 



or 



, 32\xh f 

w dw- L udx. 

E' 

Integrating over a distance Ax obtains 



(x + Ax)- 



w 



W+ J p Udx 



(6-66) 



(6-67) 



(6-68) 



If the terms under the integral can be assumed to 
be constant, this simplifies further to 



w(jc + Ax) = 



3/ X 96 ^ h f 

w (x)-\ -uAx 



(6-69) 



If the height is not constant and the fluid is non- 
Newtonian, a similar equation can be written for 
the cross-sectional area of the fracture by using the 
power law rheological parameters: 



where 



(A) + - f — u"Ax 



-il/Cn+2) 



<D = - 



1+2 1 

ef w(y) \ n. 
\\\ w J 



dy. 



, (6-70) 



(6-71) 



h f (2 + 4n)J\ 

For an analytical solution, Ax would be the entire 
fracture length (Nolte, 1991), and this would be 
combined with a tip criterion and a volume-balance 
equation. The numerical solution proceeds simi- 
larly, except that Ax is chosen sufficiently small to 
obtain an accurate solution. Fluid loss is integrated 
over the time step, which allows obtaining accept- 
able accuracy, even with large time steps. The solu- 
tion method at each time step is as follows: 

1 . Estimate a tip velocity. 

2. For each element, working in from the tip to the 
well, 

a. calculate an average fluid velocity based on 
the velocity at the outer side of the element 
and the estimated velocity at the inner side 
(At the first iteration, assume the inner fluid 
velocity is equal to the outer fluid velocity.) 
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b. determine the cross-sectional area at the inner 
side of the element for the estimated velocity 
by using Eq. 6-70 

c. determine the velocity at the inner side such 
that the leakoff and element volume change 
during the time step result in a mass balance 

d. repeat the iteration using the new velocity. 

3. Compare the actual flow into the fracture with 
the wellbore velocity calculated by the iterative 
scheme in the preceding step. 

4. Refine the estimate of the tip velocity using a 
Newton-Raphson method until volume balance 
is achieved, which typically takes two to four 
iterations. 

This method of solving the equations is efficient 
because the velocity does not vary significantly 
along the fracture for typical cases. For typical 
PKN cases with a single fluid, the fracture can 
be divided into about 10 elements. For non-PKN 
cases, the grid must be chosen sufficiently fine 
that the integrand in Eq. 6-68 (which includes 
effects of fluid rheology and fracture height) is 
approximately constant in each element (because 
the solution scheme is derived with the assump- 
tion that it is constant). 

Regardless of whether a moving-, or fixed-grid 
method is used, usually only a small number of 
elements (about 10) is necessary to obtain a 
reasonably accurate solution to the equations 
described so far. However, other information may 
be required at a much finer resolution. To achieve 
this, the schedule is typically divided into a large 
number of substages (about 100). Quantities such 
as proppant concentration, fluid temperature and 
acid concentration can then be tracked on this 
finer grid. In addition, particularly in acid frac- 
turing, it is desirable to track leakoff and etching 
on a finer grid. To do this for methods using 
a moving grid, a second grid that does not move 
is established. Quantities such as reservoir tem- 
perature, proppant bank height and leakoff vol- 
ume in the reservoir are tracked on this solid- 
based grid. 

Nonequilibrium-height solution 

It was noted in "Solid mechanics solution" in 
Section 6.3-2 that the assumption of slow height 
growth allows creating a pressure-height-width 
table prior to solving the equations of fracture evo- 
lution. This so-called equilibrium-height assump- 



tion is quite accurate, provided that the fluid is 
moving relatively slowly in the vertical direction so 
that the pressure drop resulting from vertical fluid 
flow is negligible. This assumption is violated if 
high-permeability zones are exposed, because fluid 
must then move rapidly because of the increased 
leakoff in such layers. Also, if the stress in the sur- 
rounding zones is insufficient to confine the frac- 
ture and the vertical tips extend quickly, then the 
fluid must move quickly to fill the resulting frac- 
ture. In either of these cases, the pressure gradient 
resulting from vertical fluid flow may become 
large, and the equilibrium-height assumption 
becomes invalid at these locations in the fracture. 

To remove this assumption and obtain valid 
results from a simulator, some restriction must 
be placed on height growth. For nonequilibrium- 
height growth, the pressure gradient because of 
fluid flow in the vertical direction must be approxi- 
mated, based on the rate of height growth. It is 
common to base this approximation on the KGD 
model (e.g., Settari and Cleary, 1982). In Section 
6-7 on tip effects, an analytical near-tip solution 
developed by Lenoach (1994) is discussed that pro- 
vides an expression for the net pressure of the form 



3**2(8 Mil 



1*.H 



VM£2+4 



V ' (6_72) 



where u, ip is the tip velocity and p is 2/(2 + ri). As 
previously noted, for a fracture under constant pres- 
sure, the stress intensity factor is related to the net 
pressure by 



(6-73) 



The Lenoach tip solution can be used to obtain 
an apparent fracture toughness caused by the non- 
zero tip velocity by combining Eqs. 6-72 and 6-73. 
This effect can be added to the actual rock tough- 
ness, and the sum is used in Eqs. 6-46 and 6-47 
instead of the actual rock toughness to determine 
the fracture height growth. The basic algorithm 
used to move from one pair of vertical tip positions 
to another during a time step is as follows: 
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1 . Estimate the top and bottom tip velocities for 
the cell. 

2. Calculate the new fracture tip positions at the 
end of the time step, using these velocities. 

3. Calculate the stress intensity factors from Eqs. 
6-46 and 6-47. 

4. Determine the excess stress intensity factors 
(i.e., the calculated value minus the rock 
toughness). 

5. Calculate the velocities required to generate 
these excess stress intensity factors, using Eqs. 
6-72 and 6-73. 

6. Compare the velocities with the guessed values 
and iterate until the correct velocities are found. 

One of the advantages of the equilibrium-height 
models is the speed gained by precalculating a 
table of the fracture height-pressure relation. Not 
only is this not possible for the nonequilibrium 
model, but the iterative process to determine the tip 
positions can be time consuming. The nonequilib- 
rium-height algorithm should therefore be used 
only when necessary because of the apparent rapid 
height growth indicated by the equilibrium-height 
calculation. 

Lateral coupling 

In the description of the solid mechanics solution 
provided previously, the basic assumption is that 
individual cross sections act independently (i.e., 
plane strain in the horizontal direction, or laterally 
decoupled). This is implicit in the assumption that 
the pressure and width at any point are uniquely 
related. In reality, the pressure at any point is depen- 
dent not only on the local width, but also on the 
width distribution over the entire fracture, as dis- 
cussed in Section 6-3.1 on planar 3D models. This 
lateral coupling is generally not important, unless 
the fracture wing length is less than the height. Even 
then, the fracture geometry will not be significantly 
different if lateral coupling is neglected, although 
the pressure response may be underestimated. Lat- 
eral coupling can be included in the solutions 
described previously (see Sidebar 6E). 

The effect of lateral coupling during pumping is 
to increase the pressure at and near the well and to 
decrease it near the tip. Figure 6-11 shows the evo- 
lution of pressure during a treatment for a confined 
fracture simulated using the KGD, PKN and later- 



ally coupled PKN models. The pressure predicted 
by the laterally coupled model is always higher 
than either the KGD or PKN solution would pre- 
dict. It can also be shown that the width at the well 
is always smaller than that predicted by either of 
the simple models. The point in Fig. 6-11 where 
the pressure from the laterally coupled model is 
lowest (and where the pressures from the KGD and 
PKN models are equal) corresponds to a square, 
where the fracture wing length is one-half of the 
height. The pressure calculated by the laterally cou- 
pled model exceeds that predicted by the KGD or 
PKN model at this time by approximately 40%, 
which is comparable to the pressure in a radial 
fracture of similar dimensions. 





1400 




1300 






8? 




D 

8 


1200 


Q) 




0. 






1100 




1000 



1 










i 










\ 






























— * 






Lateral coupling 

KGD 

PKN 





























10 15 20 
Time (min) 



25 



Figure 6-11, Pressure record with and without lateral 
coupling. 



6-3.3. Lumped pseudo-three- 
dimensional models 

Lumped models are an alternative to cell-based mod- 
els and were first introduced by Cleary (1980b). 
Although more details are presented in subsequent 
paragraphs, it is worthwhile at this point to quote two 
sentences from the conclusions of his paper: "the heart 
of the formulae can be extracted very simply by a 
nondimensionalization of the governing equations; the 
remainder just involves a good physico-mathematical 
choice of the undetermined coefficients" and "results 
could be presented in the usual format of design 
charts, based on dimensionless groups extracted, . . . 
[a] more appealing procedure may be to program the 
solutions for a suitable pocket calculator, with the sep- 
arately determinable y or T coefficients and job para- 
meters as input." Although numerous papers have 
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been presented on the use of this model (e.g., Cleary 
et aL t 1994), which is now on laptop computers rather 
than pocket calculators, these sentences capture the 
essence of lumped models in that they are extremely 
simple models and the key to their successful use 
is the selection of appropriate coefficients for the 
problem analyzed. 

Like all the models presented previously, the start- 
ing point of the lumped model is the set of basic equa- 
tions defining the hydraulic fracturing process, which 
are mass conservation (Eq. 6-39), the relation between 
the distribution of crack opening over the length of 
the fracture 2L, net pressure distribution (similar to 
Eq. 6-36) expressed as 

L 

pM=\'{x,x')w{x')dx\ (6-74) 

-L 

and conservation of momentum (Eq. 6-40) expressed as 

qq m - ] =-jy n ~ m+] Vp/\l, (6-75) 

where y 4 is the channel factor (Vn for a Newtonian 
fluid), and various combinations of the power law fac- 
tors m for turbulence and n enable consideration of 
both non-Newtonian fluids and turbulent flow. 

In the lumped models, these equations are simpli- 
fied by assuming a fracture shape and adopting a spa- 
tial averaging approach to reduce them to ordinary 
differential equations in time. This approach implicitly 
requires the assumption of a self-similar fracture 
shape (i.e., one that is the same as time evolves, 
except for a length scale). The shape is generally 
assumed to consist of two half-ellipses of equal lateral 
extent, but with different vertical extent. 

It is instructive to consider some of the lumped 
equations for the KGD model (Cleary, 1980b). The 
mass balance is obtained by averaging over the frac- 
ture length: ' 

p(qw-Lq L ) = d(y 2 pwL)/dt 9 (6-76) 

where 

"~^PJE)L (6-77) 

q m w m =y;(wf n + 2 /L\ (6-78) 

where 

Y^Y 2 Y 4 £/M, ( 6 " 79 ) 

which is the ID form of Eq. 6-75. Equation 6-76 is 
similar to Eq. 6-15 (based on Carter, 1957) with the 



addition of y 3 , and Eq. 6-77 is identical to that of 
Geertsma and de Klerk (1969) with y 3 replacing 
4(1 - v 2 ). Superficially, these equations are extremely 
simple, but the values of the y coefficients are not 
always obvious and may not be constant. As noted 
by Crockett et al (1989), these models are extremely 
general, with the degree of accuracy limited ultimately 
only by the effort invested in determination of the y 
coefficients by detailed simulations, laboratory experi- 
ments or field studies. 

For more general fracture shapes (i.e., with height 
growth), it is typically assumed that height growth is 
governed by a KGD-type solution and length growth 
by a PKN-type solution (Cleary et al., 1983), although 
this is not a theoretical limitation of lumped models. 

One area in which lumped models have been 
exploited extensively is in the development of com- 
puter software systems to apply and use pressure data 
during a treatment. Some of the key characteristics 
and requirements for such a software system are that 
(Crockett etal, 1989) 

• the physics is realistic and general 

• execution time is much faster than treatment time 
to allow repetitive execution for pressure history 
matching 

• the software can use improved estimates of parame- 
. ters obtained in real time (i.e., during the treatment). 

Although these real-time software systems are gen- 
erally referred to as real-time hydraulic fracture mod- 
els (e.g., Crockett et al, 1989), the model itself is only 
a small part of the software and should address the 
first requirement listed (i.e., realistic and general 
physics). The second and third requirements are com- 
puter hardware and software design constraints. 
Because lumped models were developed for pocket 
calculators in 1980, they impose a minimal impact on 
computer hardware and systems. As computing power 
continues to improve, it will become possible to run 
increasingly sophisticated models during treatment 
execution either at the wellsite or remotely. There are 
other software design issues, such as robust execution 
with a wide variety of parameter values, easy import 
and superposition of actual data on model output, and 
graphical display, that are required for a useful soft- 
ware system for real-time applications. Discussion of 
these issues is beyond the scope of this volume. 
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Analysis of a Vertical Crack in a 
Multilayered Medium 

A general mathematical formulation to analyze cracks in a multilayered medium is 
constructed. First, a matrix for a single layer is formed in the Hankel transformed 
domain. Then a global matrix is formed by assembling together each layer matrix. 
The displacement and stress fields are obtained by inverting the Hankel transform. 
Finally, the planar crack problem is solved by the boundary integral method. The 
results given are the crack opening displacement and the stress intensity factors. 



1 Introduction 

Hydraulically-induced fracture is a technique being 
developed in the petroleum industry. The basic idea in the 
hydraulic fracture process is to use a highly pressurized fluid 
to break the rock so that the resources, oil or natural gas, can 
be more easily extracted. One of the difficult analytical prob- 
lems for such fracturing processes is how to calculate the im- 
portant quantities of physical interest when the process zone 
extends from the pay zone into upper and lower barriers. The 
methodology applied in this article, which studies one layer 
having upper and lower half space barriers, can also be ap- 
plied to more general multilayered structures. 

Multilayered structures are often encountered in many 
science and engineering branches such as geophysics, 
laminated composite materials, and highway pavements. The 
first systematic approach to multilayered media may be traced 
back to the Thompson (1950) -Haskell (1953) methods. Later, 
Gilbert and Backus (1966) introduced the propagator matrix 
method; Fuchs (1968, 1971) introduced the reflectivity 
method, and Pao and Gajewski (1911) have written a long ar- 
ticle about generalized ray method. 

The Green's functions in a multilayered media were derived 
by Luco and Apsel (1983). Xu and Mai (1987) give a simpler 
formulation and show that the problem can be solved by a per- 
sonal computer. Kundu and Mai (1985) discussed the difficul- 
ty associated with the numerical integration in the trans- 
formed domain. Small and Booker (1986) solved a 
multilayered system by the flexibility matrix method. Jn Sec- 
tion 2, a general scheme to analyze the multilayered system is 
developed. 

For the three-dimensional crack analysis in an unbounded 
medium, the boundary integral method is still preferred since 
the unknowns and equations occur only on the crack surfaces. 
The unknowns in the boundary integral method can be the 
crack opening displacement, Murakami and Nemat-Nas'ser 
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(1983), the first derivatives of crack opening displacement, 
Weaver (1977), or the second derivatives of the crack opening 
displacement. The choice depends on whether or not it is more 
convenient to avoid the high order singularity. Recent articles 
by Lee and Keer (1986) and Lin and Keer (1987) use the 
triangular element in the same manner as Murakami and 
Nemat-Nasser. Although the triangular element is somewhat 
easier for mesh generation, a four node element which may be 
generated via conformal mapping, can evaluate the finite pan 
integral more easier. Other articles such as Krenk and Schmidt 
(1982), Martin and Wickmen (1983), Lin and Keer (1986), 
etc., are concerned with the scattering due to a crack in the fre- 
quency domain. 

An extra difficulty encountered for analyzing a vertical 
crack in the multilayered medium is the point singularity that 
occurs at the intersection of the crack edge and the interface. 
Several articles, Benthem (1977), Bazant and Estenssoro 
(1979), and Somaratna and Ting (.1986), have discussed this 
type of point singularity. 

In Section 3, the boundary integral equation will be con- 
structed by a Green's function approach. Triangular elements 
and a square root singularity along the crack edge will be 
adopted. Although,, the point singularity at interfaces is not 
considered in this article due to complexity, the interfacial 
phenomenon is discussed in the numerical results. 

2 Mathematical Formulation for Multilayered 
Medium 

To analyze a multilayered medium systematically, there are 
three basic principles involved. First, the matrix for' a single 
layer is formed in such a way that the same formula can be ap- 
plied to every layer. Next, the interfacial continuity conditions 
are executed in the transformed domain instead of in the 
spatial domain, since the mathematical operation required in 
the transformed domain is algebraic. The complexity of ex- 
ecuting continuity conditions in the spatial domain can be 
found in Chan et al. (1974). Finally, to reduce the computa- 
tional work, the solution is obtained by inverting the Hankel 
transform instead of inverting the double Fourier transform. 
Using these principles, the methodology to analyze the 
multilayered medium is constructed. Essentially, three types 
of matrices can be formed, namely the flexibility matrix, stiff- 
ness matrix, and propagator matrix. However, only the flex- 
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ibility matrix is considered in detail. The other two methods 
are described in Appendix A. 

First, consider a linear elastic isotropic layer which occupies 
the spatial domain 0<z<H. Although initially there is no 
body force, such forces will be introduced later into the 
system. The displacement field in this case may be expressed in 
terms of Papkovich-Neuber potentials in rectangular coor- 
dinates. However, these potentials lead, unavoidably, to a 
double Fourier transform; therefore, the following alternative 
forms in cylindrical coordinates, Green and Zerna (1954), are 
considered: 



u z = <p :i . + z<t> rtZ -K<t>, 



(Iff) 
(\b) 
(lc) 



Here k = 3-4m, v is Poisson's ratio, <p., <f> r and <j> 6 satisfy 
Laplace equations, and the comma denotes differentiation. 
These forms are particularly useful when the applied body 
force is a vertical or a horizontal point force since only the 
constant term or the first term in Fourier series expansion is 
needed. Next, the axially symmetric and asymmetric cases are 
considered. 

(A) Axially Symmetric Case. In the axially symmetrical 
case, w 5 = 00=O and <f>. and 0, have no angular dependence. 
Thus, the potentials can be solved by using the zeroth order 
Hankel transform. The results are exponential functions with 
unknown coefficients in the transformed domain. The 
displacement field is obtained by replacing the potentials in 
equation (1) with the Hankel integrals, and the stress field is 
obtained by -further application of Hooke's law. Explicitly, 
the results are: 



(2a) 



(u z> o zz )=\ (e/„,S K )f/ 0 ({r)</€ 



(2b) 



and 
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>4,exp(-{*) 
/4 2 exp(-£z) 
Aitxpti(z-H)] 
l^exptfU-//)], 



(3) 



Here, p is the shear modulus, a = 1 - 2v, and /!,(£) are func- 
tions to be determined from the boundary conditions. It 
should be noted that symbol t/„, instead of £/.,'is used in 
equation (2a); U, will be reserved for the more lengthy asym- 
metric case. 

To form the flexibility matrix, three vectors are defined as: 
(u) 7 *- ( '-U a (H) t -U r {H) y t/ a (0), t/ r (0)] (4a) 

k I T = 1 S„ (H) , (H) , 5^(0), S n (Q) ) m 

I A| r = \A^A 2 ,A } ,A A }. (4c) 



To form the stiffness matrix or propagator matrix, (u) and 
1 a\ are defined in other forms. Details about the two other 
matrix methods are given in Appendix A. After defining (u | , 
la], and (A), the (u) - (A) and \ a\ - | A) relationships are 
constructed by substituting z = 0 and z = H into equation (3) in 
the form of equations (4g)-(46). The results are written as: 



M=[f,](A) 

M = [f 2 ](A). 
The flexibility matrix for the layer is: 

[U=[f l ][f2]' 1 

[U'M = |u|. 



(5a) 
(5b) 

(6b) 



If a body force is applied in the layer, the total displacement 
and the stress fields can be written as: 

f")r=luiF+(u) (7a) 

Mr= Mf+ M- (lb) 

Here, the subscript "F' denotes that the field is due to the 
same body force as applied in the "full" space case, and 'T ' 
denotes that the field is contributed from both the body forces 
and the boundary conditions. By substituting equations (la) 
and (lb) into equation (6b), one obtains the result: 

[UMr= ("ir+[Uk)F-(uU (8) 

Next, the displacement and stress fields, [u} F and [a] F , due 
to a unit vertical point force applied at x 0 will be constructed. 
In this case, the displacement field in the spatial domain is 



Ujr = C[K/R+ ( z -z 0 ) 2 /R*] 
u fF = Cr(z-z 0 )/R\ 

Here, 

C = \/\6irp(\ -u) t R = \x-x 0 \,r=[R 2 -(z-z 0 ) 2 V /2 

The displacement and stress fields in the transformed domain 
are obtained, by performing the related Hankel transforms as 
given in equations (2a)-(2b). The results are: 

C/ aF = C[ic/f+ l*-« 0 l]exp(-{ \z-z 0 \) 

U rF = C(z-z 0 ) exp(-flz-jj) 

= -2 M C[(2-2i0sgnfe-* o ) 

+ «(*-z 0 )]exp(-f \z-z 0 \) 

S rzf= -2nC[\-2v + t \z-z 0 \) exp(-£ \z-z 0 \). 

Then [a\ F and (u]^ in equation (8) are obtained by 
substituting z = 0 and z = H into equations (\0a)-(\ \b) in the 
forms of equations (4o)-(4d). 

If the layer is replaced by an upper half space, then 
/1 3 = A 4 =0 in equation (3), and there is no change in the for- 
mulation. Hence, the half space solution may be viewed as a 
special case of a layer solution. 

For a multilayered system that contains N layers, the re- 
maining procedure is to assemble the matrix of each single 
layer to form the global matrix. To demonstrate the assembly 
process, two consecutive layers / and /+l are considered. If 
layer / is above layer /+ 1, the interfacial conditions are: 

( S^O. ) , 5^(0, ) ) r = 1 S :z (H l+l ) t S rz (H i+ , ) ) T (\2a) 
IVzzW&UAOMr^lU^Hi+^UAHi^nr. (\2b) 



(9a) 
(9b) 

(9c) 



(\0a) 
(10b) 

Ola) 

(\m 
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The flexibility matrix after assembling is: 
' [F tf ] = [fJ, + [U +l . 
The transposed vector of | a \ T is: 

I S a (//,), S n (Hi), S tt <0,), 5 r ,(0,), SJ0 UI ) 9 5, : (0 ( .,,)] t- 
The transposed vector of [ u ) T is: 
[-(/ a (J/,). -l^(// ; ),0. ) 0.,U B <0, + I ) l i/,(0 /+l )] r . 
The forcing terms are: 
(IMIal/r-lul^+dfJIalf-lul^,. 

• (B) Asymmetric Case. After discussing the process to 
solve the axially symmetric case, the displacement and the, 
stress fields due to a point ' 'horizontal' 1 force are considered. 
Similarly to equation (8), the layer equation is: 

Here, the goal is to construct the matrix [f A ], and vectors 
\a h \ F and \u h .] F in equation (13). It should be kept in mind 
that [f h \ depends on the definition of [u h ] and |<rj. 

It is well known that the displacement field due to a unit 
horizontal point force applied at x 0 in the full space is: 



= Cr(z-z 0 )cos8/R 3 
u rF = CWR + r 2 /R l ]zo$Q 

U eF = -CK$\T\d/R 



(14*) 
(146) 
-(14c) 



where r cos0 = ;c-x o , n\nd = y-y 0 . Thus, the angular 
dependence for the potentials in equations (1) are 

Wj, 4> n <t> 9 ) = (*-"cos0, $ f cos0, $ tf sin0), 

which are solved by the first order Hankel transform; the 
results are: 

jjl^expc-^j+^acxplftz-//)])^^)^ (15) 

where ($, , $ 2 , $3) = Bv substituting equation (15) 

into equations (1), the displacements field in terms of Bj is ob- 
tained. The stress field is obtained by differentiating the 
displacement field. The results are in the following compact, 
forms: 

(w,/c, w,/c + Mf,/5, u r /c-u e /s, 

o zz /c, o n /c + o 6z /s t o^/c- o 6z /s\ 



-1 



[tV,({r), U p Jj{ir), U m J 0 (Hr), 



To obtain the flexibility matrix, ( u A ) , (a A | and j B A j are 
defined as: 

\u h ) T = [-U Z [H), -U P {H), 

-U m (H), U,(0), U p (0), U m (0)) (18a) 

\o„\ T = iS : {H), S,{H), S„(H), S : (0). S p (0), S m (0)| 

(186) 

(BJ T = |B,, B 2 , B 3 , B 4 , B 5 , B 6 ). (18c) 

Then, the (uj - (B A ] and \ a h | .- |B A ) relationships are ob- 
tained by substituting z = 0 and j = W into equations (17). The 
results are written as: 



|uJ=[f 3 ]|B A ) 
where the flexibility matrix is 



(19a) 
09b) 

(20) 



SJ^r)rS p J 2 {ir) t S m J 0 ar)}^ 
where c = cos0, s = sinfl and 



(16) 



To obtain \u h ] F and \o h \ F > the appropriate Hankel 
transforms as indicated in equation (16) are performed on 
equation (14). The results are: 

C (z-z 0 ) exp(-£ lz-z 0 1) 

(V= C (1/?+ Iz-zJ) ' exp(-£lz-z 0 l) 

ty mf = C[(7-8.)/£- iz-z 0 l] exp(-£lz-z 0 D (21) 

5^/2 M - C(l-2*-£ lz-z 0 !) exp(-f lz-z 0 l) 

V-/2 M = C[-£(z-z 0 )l exp(-£iz-z 0 i) 

S mf /2 M = C[( - 4 + 4 I /)sgn(z - z 0 ) + £ ( z - z 0 )]exp( - £ lz - z 0 1 ). 

Then, |u^] r and { a A ) /r in equation (13) are obtained by 
substituting z = 0 and z = H into equation (21) in the forms of 
equations (18fl)-(18c). 

Generally speaking, when the loading of boundary condi- 
tions are more complex, more terms for the Fourier series ex- 
pansion may be needed. However, the formulation is basically 
the same as the first term expansion, and the main differences 
between the mth term and the first term are in the replacement 
of 7 0 (£r), y,(£r), J 2 (£r) by J^.tfr), J„Ur), J^Ur). 

3 Integral Equation and Numerical Scheme 

To set up the integral equation, two symbols, C ik {x, x 0 ) and 
T Uk (x, x 0 ), are introduced, where G ik (x, x 0 ) is a representa- 
tion for the /th direction displacement at x due to a unit point 
force applied along kih direction at x 0 and r ( ^(x, x 0 ) is the 
•corresponding ij stress component. The processes required to 
derived the Green's functions have been described in Section 
2. For a vertical planar crack with surface normal to x axis, the 
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Fig. 1 The mesh and the coordinate system ol the numerical examples 
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Fig. 3 The normalized crack opening, M Au 7 /p<2ae) 1/2 ( along the edge of 
the crack, angle 7 = tan " 1 (z/y) for bimaterials 
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Fig. 4 The normalized crack opening, ^Au 7 /pa, along the Z-axis for 
bimaterials 
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Fig. 2 The normalized stress intensity factor, K,/K tF , along the edge of 
the crack, angle 7 = tan " 1 {zJy) for bimaterials 



integral representation for the displacement field due to the 
mode I crack opening displacement is: 



w *(x 0 ) = W uk (x; x 0 ) &u^x)dA (x). 



(22) 



In Appendix B, it is shown that equation (22) can be derived 
from the equivalent body force method. 

The integral equation is obtained by applying Hooke's law 
to equation (22) and applying the pressure boundary condition 
along the crack surface. The result is: 

-P(x 0 ) = a n {x 0 ) = \K{x,x 0 )&u x (x)ctA{x) (23) 

where 

^(x t x 0 ) = 2 Mo ((l-r 0 )r 11Uro + , 0 [r ll2i>0 + T M ^ 0 ] )/(l-2^) 

(24) 

and n 0 = n(x 0 ), v 0 = \>{x 0 ). Details concerning K(x t x 0 ) are 
given in Appendix C. Numerically, equation (23) is written in 
discrete form as: 

P(x 0 )=f ( K(x; x 0 ;OA«,(x)w(x)£7V4(x)l (25) 

where w(x) is the weighting function for the crack opening 
displacement, i.e., Au, (x) = w(x)Au l (x), K(x, x 0 ; £) is the 
kernel in J domain, and S m is the subdomain. The numerical 
integration over the S m domain is described by Lin and Keer 
(1987). 



4 Numerical Examples and Discussions 

The mesh an d the coordinate system of the numerical ex- 
am p>s^ajiIiI}0N|n in Fig. 1. The crack is a penny-shaped 
crack, fadius cy>and the materials consist of an upper half 
space, a^fttpaflayer, and a lower half space. The central layer 
represents the pay zone, while the upper and the lower half- 
spaces represent barriers. The exact properties of the materials 
and the locations of the two interfaces are specified in the ex- 
amples. The loading is by a constant fluid pressure p y applied 
on the crack surfaces. 

In the first example,..two interfaces are. located at z = 0 and 
z = 0. 3a. The. shear moduli of the upper half space, the central 
layer and the lower half space are 4/x, 4/x, and ji, respectively. 
The Poisson's ratios are 0.2 for all materials. For this case, the 
kernel /Trias a closedform (Lee and Keer (1986)) which can be 
used to check the numerical integrations of the Hankel 
transforms. After checking, the Hankel transforms are per- 
formed by Simpson's rule with 400 integration points and with 
the maximum £ = 160/ff. In Fig. 2, the normalized stress inten- 
sity factor along the crack edge is plotted. The normal- 
ized factor K !F = 2p{a/ny n corresponds to the stress intensity 
factor for a penny-shaped crack in fully space subjected to 
uniform tension. The result is same as in a previous paper by 
Lin and Keer (1988). It is noted that the stress intensity factor 
is greater than K IF in the upper medium and is smaller than 
K IF iri the lower medium. The stress intensity factor at or near 
the interface is arguable because the point singularity at the in- 
tersection is not taken into account by the weighting function. 
In Fig. 3, the crack opening is plotted along the edge. The 
curved line gives the crack opening divided by the weighting 
function, (2ae) w2 > where e is the distance from the collocation 
point in the boundary element to the crack edge. To under- 
stand this curve, two extreme cases are compared. When the 
shear modulus in the region z>0 is replaced by then, from 
the full space solution, ^Au l /p(2ae) l/2 = 4(1 - which is an 
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Fig. 5 The normalized crack opening, fi^u 1 /p (2ae) , along the edge 
of the crack, angle -y = tan ~ 1 . (z/y) tor trimaterials 



Fig. 7 The normalized crack opening, ^JkU 7 /pa, along the /-axis for 
trimaterials 




upper bound for the crack opening (upper dashed line). When 
the shear modulus for z<0 is replaced by 4/i, then, 
jiAW|/p(2£re) l/2 = (l-t>)/7r, which is a lower bound for the" 
crack opening (lower dashed line). The crack opening along 
the edge is decreasing from the lower medium to the upper 
medium because the upper medium is suffer. Rapid changes 
occur near the interface. In Fig. 4, the crack opening along 
^ = 0 is plotted. The two dashed lines correspond to the upper 
bound, M Au l //70 = 4(l-^)[l 2 -(z/fl) 2 ] i/2 /7r, and the lower 
bound, pAu^/pa^ (l - v)[\ 2 - {z/a) 2 ] W2 /ir, of the crack open- 
ing as in the previous discussion. 

In the second example, two interfaces are located at z-0.3a 
and z- -0.3<3. The shear modulus is p for the central layer 
and rj/i for the upper and lower half spaces. The Poisson's 
ratio is 0.2 for all materials. The three curves shown in Fig. 5 
are the opening around the edge divided by the weighting 
function for jj = 2, 4, and 10. The lower bound for these curves 
are pAu x /p{2ae) 1/2 = 4(1 - v)/t)-k = 1 .02/17. The three curves in 
Fig. 6 are the corresponding stress intensity factors. The con- 
clusions are the same as the those obtained in the previous ex- 
ample: the opening in the softer medium is larger, rapid 
changes in crack opening occur around interfaces, and the suf- 
fer region has a higher stress intensity factor. Figures 7 and 8 
are the corresponding crack opening along Y axis and Z axis, 
respectively. 

Remarks: First, if there are many layers in the system, the 
method of propagator matrices given in Appendix A is recom- 
mended. The flexibility matrix is used here only because there 
are only two interfaces. Second, the current approach can be 
applied in the frequency domain (see Xu and Mai (1987)). 



Even for the transversely isotropic medium, the current ap- 
proach is also applicable (see Small (1986)). 
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APPENDIX A 

Stiffness Matrix and Propagator Matrix 

Al Stiffness Matrix: To obtain the layer stiffness 
matrix, |u) and \<x\ are defined as: 

\u} T =\U zz (H), U,(H) t UJ0), U r (0)\ (Ala) 

l*) T ={-S s {H) t -S n {H), S a (0), 5^(0) ] . (Alb) 

Then the ( u | - ( A ] and ( a \ - ( A ) relationships are obtained by 
substituting <; = 0 and z~H into equations (3). The results are 
written in the same forms as equations (5a) and (5b)\ the stiff- 
ness matrix is: 

M^Kf.]- 1 - (A2) 

If a body force is considered, the equation for the layer is: 

[sj|u| r = k) 7 +[sJ(u) F - \a\ F . (Al) 

Where \w) F and \u] F can be obtained from equations 
(10fl)-(l Wand (A\a)-(A\b) by substituting z = 0 and z = H. 

For the case of applying a horizontal point force, the ( u h ) 
and ( o h I are defined as: 

(uj r = i U z (H) i U p (H), U m [H), U z (0) t t/ p (0), U m (0)\ 

(A 4a) 

\o- h \ T =\ -S.(//). -S p {H) t -S m (H) t S z (0),S p (0),S m (Q)) 

(A 4b) 

Then, the corresponding [sj, [o h \ F and ( u h } F can be found 
by applying equations (17), (21) and (A4a)-(A4b). 

A2 Propagator Matrix: To form the propagator matrix, 



two vectors |T]and |B) for the top and bottom surfaces of a 
layer are defined as follows: 

( T ) T = [ U a (H) , U, (H) , S a (H) , S^H) \ (A 5a) 
I B 1 7 = { U u (0), U, (0), S a (0), S rz (0) 1 . (A5b) 

By substituting z = 0 and z-H into equation (3) to find the 
(T)-(A) and (B|-(A| relationship, the results are: 

|Tl = [ gl ][A] (A6a) 

|B) = [g 2 ][A). (A6b) 

The propagator matrix [P] is the relationship between (T ) and 
IB), thus, 

[P][T) = (B) (Ala) 

[P] = [g2][g 1 ]- 1 . (Alb) 

To demonstrate how the propagator matrix works, two con- 
secutive layers, / and /'+ 1 are considered. Then, 

[P]/(T),= (B| ( . (A%a) 

[P]/.iiT) / + I = [B),. + 1 . (AZb) 

These two equations may be considered into one equation by 
applying the interfacial continuity condition, 

(B),= (T] f+l . (A9) 

The result is 

[P]/ + i[P]/[T) y =(B) ( . + l . (>t 10) 

Thus, if the system contains n layers, the relation between 
the top and the bottom surfaces are constructed as: 

[PUystemMm^lB),. (AW a) 

[Ptsystem^lPUP],,., [P] 2 [P]|. (AMb) 

Now, applying boundary conditions on the top surface of the 
first layer and on the bottom surface of the last layer, the re- 
maining unknowns in (T) l and {B)„ can be solved. After 
j T ) , is solved, the remaining ( T j are obtained by 

.iT|,= (B|,. l =[P] ( _ 1 (T] ( _ 1 . 

In the above process, the [P] matrix is the downward prop- 
agator matrix. Alternatively, the upward propagator matrix in 
a layer also can be found by defining: 

(PJ(B) = [T) (A\2a) 

|PJ = [g l ][g2]" 1 . (A\2b) 
Then, the system upward propagation matrix is 

[P u system] IB],, = [T) t (Alia) 

[P„ system] = [PJ,[P U ] 2 [PJ n -, [P w ]„. (-413*) 

Now, consider the case that a body force is present in the y'th 
layer, the equation for this layer is: 

IV\j\T)tj=1*\tj + 1*]j\T) F j-IB\ F j. (A\4) 

The vectors (T]^ and IB] FJ are easily obtained from equa- 
tions (10o)-(i \b). Next, both the downward and upward prop- 
agator matrices are applied to construct the relationships be- 
tween the top and bottom surfaces. The result is: 

(p]>[p),-, . . • [Pi.m^tpj,.,, .... [pj n iB) n 

+ [P] ; .[T)^.-(B] F; , (.415) 
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For the asymmetric case, the upward and downward prop- 
agator matrix are formed. in the same manner. (T] and (B] 
are defiheci as: 

Ml 6a) 

( B J r = ( t/ : (0), U p (0), U m (0), S : (0), S p (0), S m (0) ) 

(A\6b) 

From equations (17), (21), (A\6a-b), the propagator matrices 
can be obtained. 



APPENDIX B 

Derivation of Integral Representation From Equivalent 
Body Force Method 

When there is unbalanced strain, t mn) on a system, the 
equivalent body force on this system is ~C pqmn e mnqi where 
c pqmn are elastic constants. This strain e mn is also called 
eigenstrain by Mura (1982). Hence, the displacement field due 
to this body force is: 

M*o) = "\C pqmn e mniQ (x)G kp (x 0i x)dV(x) (B\) 

where G kp (x 0 ,x), (Section 3), is the displacement along the k 
^direction at x 0 due to a p direction unit force at x. Because 
G* P (x 0 ,x) = G M (x,x o ) due to Betti's reciprocal theorem and 
Cpqmn = C mnpQ by the symmetry properties of the elastic con- 
stants, equation (£1) is written as: 

u k (x 0 ) = \e mntQ {x)C mnpq G pk <x,x 0 )dK(x). (B2) 
Integrating by parts, the result is: 
M* 0 ) = -\tmn{*)C mnpq G pkq {x,x 0 )dV{x) 

= \t mn {x)r mnk {x,x 0 )dV{x). (£3) 

The meaning of i mnk is also stated in Section 3. Finally, e mn is 
integrated along the normal of the crack surface to obtain the 
crack opening; the result is the integral representation in equa- 
tion (22). 



APPENDIX C 

Kernel 

To evaluate /f(x,x 0 ) in equation (24), first consider r, ' 
Since x = x o = 0, 11Uo ' 

r ii3 = r m =2 M [(1 - v)u r /r+ v{ u r r + u z . )]/(l - 2u). (CI) 
Substituting u r and u., equations (2<z)-(3), the result is: ■ 
T « 3 = \o 12m? 2 /(l - 2u)] 1 2*H 2 exp( - 

+ A^expttz - $//)]/ a ( (r) + L Z J X ({r)/{ r ] tf{ (C2) 

and 

I 3 = [->!, + (l-{zM 2 ]exp(-{*) 

+ I^3 + (l+^)^4]exp(£2-f/0. (C3) 

To evaluate t W3iS0 , first evaluate (A),. 0 . This evaluation is 
achieved by differentiating ( a ) F and ( u f F with respect to z 0 in 
advance, see equations (4), (\0a)-(\\b). 

The t mi and r ll2 terms are a bit more difficult. However, 
because the layers are horizontal and material is isotropic 
within each layer, r 1Ui , 0 = - r !Ut x% r M2i v0 = - r u2 v and 
T n2(0»> ; -^o.*) = r U] [y-y 0l 0, z). Next, t im and r 221 are ex- 
pressed in terms of t„, , t^, 7 M1 . Then, performing the proper 
differentiation as indicated by equations (16)-(17), the results 
are: 

+ E.[-J 0 (Zr) + J [ (tr)/tr]}dt (C4) 

p Oo 

T m.» = j o 2/x^|-3(/ p y 2 (£r)/(?r) 2 -f (25 3 exp(-^) 

+ 2£ 6 exp( - + { *) - £ z )J l (f r ) /{ r 1 rf{ (C5) 

where t/ p is given in equation (17) and 

£, = 2*[£ 2 exp( - Zz) - B s exp( - { *)]. (C6) 

When x and x 0 are in the same layer, the kernel K{x,x 0 ) con- 
tains an extra term, the full space term, i.e., 

^(x t x 0 )=^/4x(l-^ 3 . " (C7) 



Journal of Applied Mechanics 



MARCH 1989, Vol. 56/69 




SPE 80935 



Modeling of Hydraulic Fracture Height Containment in Laminated Sand and Shale 
Sequences 

Jennifer L. Miskimins, SPE, Colorado School of Mines, and Robert D. Barree, SPE, Barree & Associates, Inc. 



Copyright 2003, Society of Petroleum Engineers Inc. 

This paper was prepared for presentation at the SPE Production and Operations Symposium 
held in Oklahoma City, Oklahoma. U.SA, 22-25 March 2003. 

This paper was selected for presentation by an SPE Program Committee following review of 
information contained in an abstract submitted by the authors). Contents of the paper, as 
presented, have not been reviewed by the Society of Petroleum Engineers and are subject to 
correction by the authors). The material, as presented, does not necessarily reflect any 
position of the Society of Petroleum Engineers, its officers , or members. Papers presented at 
SPE meetings are subject to publication review by Editorial Committees of the Society of 
Petroleum Engineers. Electronic reproduction, distribution, or storage of any part of this paper 
for commercial purposes without the written consent of the Society of Petroleum Engineers is 
prohibited. Permission to reproduce in print is restricted to an abstract of not more than 300 
words; illustrations may not be copied. The abstract must contain conspicuous 
acknowledgment of where and by whom the paper was presented. Write Librarian, SPE. P.O. 
Box 833836, Richardson, TX 75083-3836 U.SA., fax 01-972-952-9435. 



Abstract 

The results of hydraulic fracture modeling in a hydrocarbon 
reservoir consisting of laminated sand and shale sequences are 
presented. Work previously published by other authors 
addresses hydraulic fracture growth characteristics in layered 
intervals, but little published information is available for 
growth in thin-bedded reservoirs. For the purposes of this 
paper, thin beds are considered to be on the order of less than 
6 inches (15 cm) in thickness. In certain reservoirs, there may 
be thousands of these thin beds involved in influencing 
hydraulic fracture growth. 

The hydraulic fracture modeling is based on data from the 
North LaBarge Shallow Unit (NLBSU) located in the Green 
River basin of Wyoming. The NLBSU produces hydrocarbons 
from the Cretaceous-age Mesaverde interval. This reservoir is 
comprised of laminated sand and shale sequences. Hydraulic 
fracturing is necessary for economic production. Actual field 
data from two NLBSU wells is used in the analysis. The data 
includes stimulation treatment pressures and information, as 
well as a radioactively tagged fracture treatment. 

The ^hydraulic fracture modeling software GOHFER v is 
-;used in;the analysis. GOHFER is used based on its ability to 
approximate decoupled rock deformation with shear-slip, ; and 
its extensive use in the industry. 

The results of this study strongly indicate that the sand and 
shale laminations present in the Mesaverde reservoir act as a 
hydraulic fracture height containment mechanism. Two events 
are presented as the reasons for this observed containment. 
First, shear slippage, is prevalent due to the numerous 
interfaces present in the subject depositional environment. 
These interfaces inhibit fracture growth through shear failure. 
Second, the contrast in rock mechanical property values 
between the sand and shale intervals reduce the fracturing 
energy available for fracture propagation. 



Further research in this area is strongly recommended due 
to the significant oil and gas reserves associated with this type 
of reservoirs. These additional research areas include 
laboratory investigations and field work with tiltmeter and 
microseismic analysis. Meanwhile, fracture modelers are 
strongly encouraged to evaluate the effects of laminations in 
their modeling efforts. 

Introduction 

Predicting hydraulic fracture growth patterns is the main goal 
of hydraulic fracture modeling. The modeling practitioner 
strives to create an input file that represents the nature of the 
reservoir as closely as possible. This paper describes modeling 
results from a depositional environment which requires a 
particularly detailed input model. The subject environment is 
laminated sand and shale sequences. 

In 1977, Daneshy 1 provided one of the first investigations 
into hydraulic fracture growth in layered formations. He 
concluded that there was a significant relationship between the 
strength of the interface bond between two formations and 
whether or not a hydraulic fracture would cross the interface. 
Although Daneshy' s work addressed layered formations, the 
results were limited to relatively thick layers of two or three 
rock types. He did conclude that slippage at the layer 
interfaces would be a relevant containment mechanism and 
that the thinner the layers, the better the chance for slippage to 
occur. 

During the same time period as Daneshy's work, \ 
Simonson, et al, published results showing the significance of 
elastic properties, in-situ stresses, and pressure gradients on 
hydraulic fracture containment. 2 Since the work of Daneshy \ 
and Simonson, et al, several authors have published treatises; 
addressing various containment mechanisms. In 1980, : 
Hanson, et al, concluded that shear stress at an interface \ 
determines whether or not a fracture traverses the interface. 3 ) 
Warpinski, et al; showed that minimal in-situ stress contrasts 
could restrict fracture growth. 4 Warpinski and Teufel 
discussed how various geological features, such as joints, , 
faults, and bedding planes, could affect hydraulic fracture'! 
growth. 5 More ;; .fecently, hv 1598;,' Barree and Winterfeld-- 
$hqwed ^how? ■ interfacial slippage- has a large influence on- 
^hydraulic ,;fracture containment through .shear decoupling. 6 ; 
Finally^ in 2001, Smith, et al, discussed how modulus in 
layered situations can be incorrectly averaged, thus providing v 
questionable results in modeling attempts. 7 
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As discussed, References 1-7 provide a brief overview of 
the main mechanisms generally accepted as containment 
mechanisms in hydraulic fracture height growth. References 
8-14 provide further considerations for hydraulic fracture 
height containment mechanisms and layered formation 
characteristics. However, there is little published information 
for hydraulic fracture growth in laminated rock formations. 
Although laminated formations are "layered," they are unique 
since there may be thousands of layers in a given zone of 
interest. These numerous layers are interactive and present 
unique problems in hydraulic fracture growth patterns. 

This paper provides hydraulic fracture modeling results for 
a reservoir consisting of laminated sand and shale sequences. 
These sand and shale laminations range from 0.5 to 6 inches 
(1.3-15 cm) in thickness with an overall reservoir thickness of 
approximately 150 feet (46 m). Detailed model develop is 
discussed along with the modeling results. Conclusions based 
on the modeling results are presented. 

Model Development 

The North LaBarge Shallow Unit (NLBSU) is the field chosen 
for use in this study. The NLBSU is located in the Green River 
basin of Wyoming (Figure 1). The NLBSU produces oil and 
gas from the Cretaceous-age Mesaverde formation at an 
average depth of 2000 feet (610 m). A variety of research on 
the NLBSU geological setting and stratigraphy has been 
performed on the area during the last several years. 15 ' 19 
Hydraulic fracturing of the Mesaverde zone is necessary for 
economic production in the NLBSU. 

Field Stratigraphy. The NLBSU Mesaverde formation 
consists of eight main subdivisions labeled C through H. A 
type log is shown in Figure 2. The Mesaverde produces from 
subdivisions F through H. Figure 3 shows an eight inch (20 
cm) section of core which displays the thin-bedded nature of 
the formation. These laminations are present throughout the 
entire Mesaverde interval and are broken down into distinctive 
sand and shale layers. Figure 4 is an FMI log from NLBSU 
#9-16E that demonstrates the laminations seen in a 270-foot 
(82 m) interval of the Mesaverde. Figure 5 is a 26-foot (8 m) 
section of this same log shown to further demonstrate the thin- 
bedded environment. 

The thin beds and main subdivisions of the Mesaverde are 
the result of a series of transgressive and regressive deposits 
generated by eustatic sea-level changes, sediment influxes, 
and crustal movement during the Sevier Orogeny. 15,20 

Modeling Software. GOHFER (Grid-Oriented Hydraulic 
Fracture Extension Replicator) is the modeling software 
chosen for this project. 21 GOHFER is a 3D finite difference 
model. Reference 22 provides a comparison of GOHFER and 
other modeling software packages. 

GOHFER was chosen for this project for a variety of 
reasons. First, the software is well known and used extensively 
in the petroleum industry for hydraulic fracture simulation. 
Second, GOHFER is a robust simulator and does not have a 
limit on the number of layers that can be incorporated in the 
input file. This was a critical consideration for the thin-bedded 
model used in this study. 



Finally, GOHFER was selected for its ability to model 
decoupled rock deformations with shear-slip. Rocks rarely 
exhibit elastic behavior, but the equations used to drive most 
fracture simulators are based on the assumption of elastically 
coupled materials. Shear slippage is known to occur during 
fracturing operations based on the occurrence of microseismic 
activity. 2 In a sense, GOHFER can allow this slippage to 
occur. Since shear slippage is believed to be a significant 
containment mechanism, 1,6 the decoupling ability provides 
significant advantages in modeling the subject depositional 
environment. 

In the GOHFER software, "process zone stress" controls 
allow the modeler to account for non-elastic behaviors, fluid 
lag, and tensile stress in the tip areas. A second input 
parameter, the "shear-decoupling radius" accounts for strain- 
energy losses at shear planes at a high angle relative to the 
direction of induced displacement (fracture width). The shear- 
decoupling reduces the full surface integral displacement 
solution so that only stresses acting over a local area affect the 
local displacement. This model drastically reduces the amount 
of strain energy concentrated at the tip of the fracture, 
removing the tip singularity and predicting improved 
containment. 

Input Models. Two wells from the NLBSU, NLBSU #9-1 6E 
and NLBSU #64-1 6ED, are used in this study. Their selection 
is based on the extensive data set available for each well. Both 
wells were hydraulically fractured during completion and 
actual field data results are available. Additionally, the 
NLBSU #64-16ED treatment was radioactively tagged and 
logged.before and after the treatment (three components of the 
fracture treatment were tagged including the pad fluid, the 
main 20/40 sand proppant, and the tail 16/30 sand proppant). 

Three modeling cases are developed for each well. The 
first case for each well is a generic base case where detail is 
not stressed. These generic cases are labeled Case 1 and Case 
2 for NLBSU #64-16ED and NLBSU #9-16E, respectively. 
Case 1 contains eight layers (4 sand and 4 shale), and Case 2 
contains eleven layers (6 sand and 5 shale). A gamma ray 
cutoff of 60% shale established the breakdown of the sand and 
shale layers. These two cases mimic the failure to recognize 
the thin-bedded nature of the reservoir. The layers in Case 1 
and 2 resemble the character of the type curve in Figure 2. 

Cases 3 and 4 (NLBSU #64-1 6ED and NLBSU #9-16E, 
respectively) are developed as intermediate cases with more 
detail than Cases 1 and 2. The individual node size for Cases 3 
and 4 are set with five- foot thicknesses, providing 
approximately 70 layers in the model. Extremely detailed rock 
mechanical property logs are available for these two wells 
from an earlier study, 17 " 18 and the GOHFER program read the 
rock property values (Young's Modulus, E, Poisson's Ratio, v) 
from a log data file based on a 1/10 inch resolution. The log 
data files are developed from the FMI log, thus providing the 
considerable resolution. 

Cases 5 and 6 (NLBSU #64-16ED and NLBSU #9-16E, 
respectively) are developed in the same manner as Cases 3 and 
4 with increased levels of detail (three-foot node thicknesses, 
providing approximately 120 layers in the model). All six 
cases are shown in Table 1. 
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With the varying degrees of detail established in the six 
cases, additional reservoir inputs such as permeabilities, leak- 
off coefficients, etc, are entered. The current fracturing 
procedure being used in the NLBSU provided the model 
inputs for fluid and proppant properties. The fluids consist of 
25 lb linear and gelled water, and the proppants consist of 
20/40 sand and 1 6/30 sand. 

Modeling Results 

The first step of the modeling procedure consisted of 
calibrating model results to the post-frac radioactive tracer log 
data from the NLBSU #64-16ED hydraulic fracture treatment. 
The authors recognize that results from tagged treatments have 
some limitations, especially if the fracture is not oriented 
along the wellbore axis, but the information provided by the 
post-treatment spectral gamma ray provides a starting point 
for calibration. Figure 6 shows further reason for confidence 
in the information provided by the spectral gamma ray. Figure 
6 shows the NLBSU #9-16E FMI log compared to the 
NLBSU #64-16ED post- treatment gamma ray. Although the 
depths of the wells are slightly different, there is a strong 
correlation between the major sand intervals of the FMI and 
the treatment placement indicated by the spectral gamma ray. 
, Significant containment of the fracture treatment is suggested 
by the spectral gamma ray. 

Using the pumping schedule, pressure information, and 
fracture height data from the actual NLBSU #64-16ED 
treatment, numerous model runs led to a modeling match for 
Case 1 for the critical data. However, this match could not be 
obtained without the addition of large process zone stresses in 
the highly laminated areas of the reservoir. The need for the 
large process zone stresses indicates significant tip effects are 
occurring during fracture growth. Additional runs led to 
similar modeling matches for Cases 2-6. During the modeling 
process, it was much easier to obtain data matches for the 
more detailed cases, Cases 5 and 6. 

Modeling of the apparent fracture height containment 
could only be obtained with the addition of the process zone 
stresses to mimic tip effects. In order to compare what the 
models would look like if the modeler did not realize or 
suspect such containment, these stresses were eliminated and 
additional runs generated. Table 2 shows the results of the 
models with and without process zone stresses. 

Leakoff was incorporated in further runs for Cases 1-6 for 
sensitivity analysis. Minimal differences observed between the 
leakoff and non-leakoff cases are most likely due to the low 
permeability of the Mesaverde formation. Efficiencies around 
50% occurred for all cases. 

Discussion 

Figures 7 and 8 show comparisons of predicted fracture 
height growth for Cases 1, 3, and 5 (NLBSU #64-16ED) and 
Cases 2, 4, and 6 (NLBSU #9-16E), respectively. Figures 9 
and 10 show comparisons of the model-predicted fracture 
half-length growth for Cases 1, 3, and 5 (NLBSU #64-16ED) 
and Cases 2, 4, and 6 (NLBSU #9-16E), respectively. 

Figures 7 and 8 demonstrate the differences in height 
containment when significantly more detail is added to the 
model input. The fracture height growth in both wells is 
reduced by 25-40% between the general (Cases 1 and 2) and 



the intermediate cases (Cases 3 and 4). The addition of layers 
to the intermediate cases had a distinct effect on the fracture 
growth pattern. Further addition of detail in Cases 5 and 6 
provides an additional effect, but it is not as dramatic as the 
change between the first two types of circumstances. When 
compared to the information provided by the tracer log, the 
input models with more detail (Cases 3-6) appear to be more 
representative of what is actually happening in the reservoir. 
The height growth of 150-175 feet suggests the fractures are 
being contained in the gross reservoir sand intervals. 

Figures 9 and 10 show the growth patterns of the model 
fracture half-lengths. As can be expected, the half-length 
growths tie closely with the height growths. For both Cases 1 
and 2 the half-lengths reach a maximum point approximately 
20 minutes into the treatment and then unrestricted height 
growth occurs. 

The modeling performed during this study and actual field 
results strongly suggest that laminations present in the 
Mesaverde reservoir are affecting fracture growth patterns. 
The main influence appears to be in fracture height 
containment. Two mechanisms provide the most likely 
reasons for the observed containment. First, the thin-bedded 
nature of the Mesaverde provides numerous interfacial planes 
along which shear slippage can occur. As Daneshy 1 suggested, 
thinner rock layers increase the chance for slippage at the 
boundaries. Also, thinner layers lead to more discontinuities 
and increase chance of slip. Additionally, the Mesaverde is a 
naturally fractured formation, which provides additional areas 
for slippage to occur. Warpinski, et al, showed the Mesaverde 
formation is susceptible to shear slippage during pumping, as 
evidenced by microseismic events. 23 

The second suspected containment mechanism is the 
contrast in the rock mechanical properties of the sand and 
shale intervals. The Young's Modulus, E, values for the sand 
and shale in the subject reservoir differ by almost a factor of 
two. This has been confirmed with both core and sonic log 
measurements. Consider now that in general if a rock exhibits 
a high Young's Modulus, the net treating pressure will be 
greater and the width will be narrower than in a rock with a 
low Young's Modulus. Equations 1 and 2 demonstrate these 
concepts for a radial model: 24 

r f 

v „" -J— — !— 

woe( *S«L ) (2. + 2) ( 2-. )( 2«2) ( 2 ) 

E 

In these equations p n is the net pressure, n is the power-law 
exponent for the fracturing fluid, K n is the coefficient for the 
power-law fluid, r f is the fracture radius, q is the flow rate in 
the fracture, and w is the average fracture width at the end of 
pumping. 

As can be seen in Equation 1, as the Young's Modulus, E, 
value increases, the net pressure will increase. Conversely, in 
Equation 2, as the Young's Modulus increases, the width will 
decrease. For the NLBSU Mesaverde, as the fracture 
propagates across a thin bed with a low modulus, the net 
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pressure decreases and the width increases. Both of these 
factors; alone and in combination, cause energy to be lost and 
fracture growth to be hindered. This is not a new idea, but in a 
reservoir with hundreds, possibly thousands, of bed contrasts, 
the energy loss can be significant. 

A related cause of shear-slip and bed decoupling is the 
stress concentration that develops at an interface between 
layers of different modulus. By Hooke's Law the deformation 
(strain) is directly related to the product of applied stress and 
modulus (Eo). If two dissimilar layers are coupled together 
and subjected to the same strain, caused by fracture width- 
opening, substantially different stresses will results in the two 
layers with a very high stress in the "stiff' layer and a lower 
stress in the "soft" material. The bonded interface between the 
layers becomes a focus of shear caused by the large 
differential stress across the boundary. In this case the shear 
slip or decoupling at the interface is caused by the induced 
deformation and does not occur except during fracture 
propagation. Once any slip occurs the tensile stress 
concentration at the tip of the fracture that allows the fracture 
to propagate as a "wedge" dissipates. 

Fracture initiation across a shear-slip interface requires that 
local fluid pressure invade the rock layer across the shear 
plane, penetrate the rock, and overcome its current in -situ 
stress and intact rock strength. This is analogous to trying to 
split a board by applying fluid pressure to its flat surface. If 
the board surface is sealed against fluid penetration there is no 
lateral surface for the fluid pressure to act on and fracture 
initiation is difficult. If the board surface can be invaded by 
the fluid then the internal "pore" pressure can be easily 
increased and the wood can be split. The fundamental change 
from mechanical or "wedge" propagation to local pressure- 
parting is a potentially large part of the fracture containment 
mechanism at sheared interfaces. 

Another contribution to shear stress development at a 
dissimilar material interface is the contribution of poro- 
elasticity through Biot's constant (a). In permeable sand the 
pore pressure is able to equilibrate during induced 
deformation. In impermeable shale the pore pressure is 
trapped and rises dramatically during deformation. This leads 
to a large drop in net normal force and helps to promote shear. 
This mechanism is again a result of induced deformation and 
does not show up as an identifiable stress or rock strength 
variation. Future work will focus on expanding the role of 
poroelasticity on fracture containment. 

The combination of shear slippage and the corresponding 
loss of strain energy within a sequence of thin beds caused by 
rock mechanical property differences and induced stress 
variations suggest that significant height containment can be 
expected in the NLBSU Mesaverde interval. But these 
characteristics are not unique to this reservoir. Any formation 
that is comprised of laminations similar to the Mesaverde has 
the potential for significant containment that is not identifiable 
through routine stress profiling. The results of this study are 
not entirely conclusive for fracture growth in thin-bedded 
reservoirs, and the authors believe additional work should be 
performed in this area. Additional mapping techniques, such 
as tiltmeters and microseismic records, could help to confirm 
the suspected containment mechanisms. Additional laboratory 
testing would also be beneficial. In the meantime, operators 



should be aware of the potential growth restrictions in thin 
beds and design their hydraulic fracture treatments 
accordingly. 

Conclusions 

Hydraulic fracture growth in laminated sand and shale 
sequences has been examined using real data from the NLBSU 
Mesaverde reservoir and the hydraulic fracture modeling 
software GOHFER. 

The conclusions from this study are as follow: 

1. Hydraulic fracture modeling indicates laminated beds 
will affect fracture growth patterns. The effects are 
manifested in fracture height containment. 

2. During modeling, significant process zone stresses 
had to be incorporated to mimic the non-elastic 
behavior of the Mesaverde formation. These process 
zone stresses account for non-elastic behaviors, fluid 
lag, and tensile stress in the fracture tip areas. 

3. A combination of mechanisms appears to account for 
fracture height containment in the thin beds of the 
Mesaverde formation. These include shear slippage, 
rock mechanical property contrasts, poroelasticity 
(drained versus undrained deformation behavior), and 
fracture blunting. Each mechanism alone may affect 
fracture propagation, but the combined effect appears 
to be significant in laminated reservoirs. 

4. Modelers should be aware of and account for the 
potential height containment caused by multiple layer 
discontinuities when designing fracture treatments for 
laminated reservoirs. If the potential containment is 
not accounted for, the entire reservoir interval may 
not be stimulated. 

5. Additional work should be done to further the 
industry's knowledge of hydraulic fracture growth, 
patterns in thin-bedded reservoirs. Significant 
reservoirs are produced from this type of depositional 
environment and further research would benefit 
reserve recoveries. This work may include laboratory 
analysis and fracture mapping techniques such as 
tiltmeters and microseismic surveys. 
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Nomenclature 

E = Young's Modulus, psi 

FMI = Formation Micro Imager 

GOHFER = Grid-Oriented Hydraulic Fracture 

Extension Replicator 



K n = 


Coefficient for the power-law fluid 


n = 


Power-law exponent for the fracturing fluid 


Pn = 


Net pressure 


q = 


Flow rate in the fracture 


rf = 


Fracture radius 


w = 


Average fracture width, end of pumping 


a - 


Biot's poroelastic constant 


v = 


Poisson's Ratio 
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in X 2.540 E+00 = cm 
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1 


NLBSU#64-16ED 


General 


8 


2 


NLBSU#9-16E 


General 


11 


3 


NLBSU#64-16ED 


Intermediate 


-70 


4 


NLBSU#9-16E 


Intermediate 


-70 


5 


NLBSU#64-16ED 


Detailed 


-120 


6 


NLBSU #9-16E 


Detailed 


-120 
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.E 2: MODELING RESULTS WITH AND 
fHOUT PROCESS ZONE STRESSES 


Case 
Number 


Fracture Half- 
Length-to-Height 
Ratio With Process 

Zone Stresses 


Fracture Half-Length- 

to-Height Ratio 
Without Process Zone 
Stresses 


1 


2.23 


0.54 


2 


4.48 


0.61 


3 


5.03 


0.64 


4 


1.75 


0.59 


5 


3.34 


0.61 


6 


2.85 


0.62 
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Figure 1: Green River basin location in Wyoming. The NLBSU is 
located at the northern extension of the Moxa Arch. 25 





Figure 3: 8-inch core section from NLBSU #9-16E. Note the 
laminated, thin-bedded sand (light gray) and shale (dark gray) 
layers. The numbers represent inch measurements and the dots 
represent half-inch measurements. 



Figure 2: Type log from NLBSU #11-17E for the producing 
Mesaverde interval of North LaBarge field. 15 
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Figure 4: NLBSU #9-16E FMI log for the depths of 1815 to 2100 ft 
(584-640 m). Note the thin-bedded nature of the formation. 
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Figure 6: Comparison of the NLBSU #9-1 6E FMI image log and the NLBSU #64-1 6ED post-treatment gamma ray. Note the containment shown 
by the tracer for the major sand zones between 1840-1990 ft. Fracture height containment is suggested by the heavily laminated upper and 
lower shale zones. 
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Figure 7: Comparison of height growth patterns for Cases 1, 3, and 5 (NLBSU #64-1 6ED). 
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Figure 8: Comparison of height growth patterns for Cases 2, 4, and 6 (NLBSU #9-1 6E). 
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Figure 9: Comparison of fracture half-length growth patterns for Cases 1, 3, and 5 (NLBSU #64-1 6ED). 



800- 















— 


































a 






W 
















V 





























































700- 



600- 



500- 



400- 



300- 



200- 



100- 



00:00 

9/26/2001 



TG Version G3.2.1 
17.Dec-02 16:45 



00:05 00:10 



00:15 



0020 
Time 



0025 



0030 0035 00:40 

9/26/2001 



Figure 10: Comparison of fracture half-length growth patterns for Cases 2, 4, and 6 (NLBSU #9-16E). 



